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CHAPTER  1 
INTRODUCTION 

1 . 1  Prel iminary  Remarks . 

The  present  global  energy  situation  has  greatly  increased  the 
demand  for  the  creation  of  more  highly  efficient  energy  conversion 
methods.  Because  of  the  preponderant  use  of  reciprocating  internal  com¬ 
bustion  engines  in  converting  the  world's  fossil  fuels  into  useful  work, 
the  need  for  innovations  in  internal  combustion  engine  design  is 
paramount . 

Until  very  recently,  the  means  by  which  such  innovations  could  be 
developed  was  limited  largely  to  experimentation.  The  theoretical 
understanding  of  the  internal  combustion  process  was  based  on  classical 
thermodynamics  and  combustion  theory.  The  properties  of  the  working 
substance  within  the  engine  were  viewed  as  being  homogeneous  throughout, 
and  the  fluid  dynamics  of  the  gases  within  the  combustion  chamber  was 
not  considered.  With  the  advent  of  great  progress  in  the  area  of  tur¬ 
bulent  modeling,  the  development  of  several  powerful  numerical  methods 
for  solving  fluid  and  energy  flow  within  enclosed  boundaries,  and  the 
significant  advances  made  in  combustion  modeling,  it  is  now  possible  to 
begin  to  analyze  with  much  more  realism,  the  processes  taking  place 
inside  internal  combustion  (IC)  engines.  The  objective  of  the  research 
effort  recorded  in  this  report,  was  to  bring  to  bear  a  select  array  of 
these  recent  developments  on  a  specific  IC  engine  innovation  in  order  to 
create  a  flexible  solution  code  for  use  in  assisting  future  experimental 
work  on  a  new  and  unique  engine  concept. 

The  diesel  engine  for  which  this  model  has  been  developed,  has  two 
significant  innovations.  A  barrier  placed  a  distance  below  the  surface 


i::-: 


of  the  piston  retards  heat  flow  through  the  bottom  of  the  piston.  As  a 
result,  the  metal  between  the  surface  of  the  piston  and  this  barrier 
stores  and  re-adds  energy  to  the  working  substance  during  the  engine 
cycle.  The  second  innovation  involves  the  use  of  a  new  combustion 
method  called  hypergolic  combustion.  In  this  method  the  fuel  is  pre¬ 
heated  to  above  its  activation  temperature,  so  that  as  soon  as  a  fuel 
radical  comes  in  contact  with  a  radical  of  oxygen  after  the  fuel  has 
been  sprayed  into  the  engine,  the  complete  combustion  of  the  fuel  radical 
takes  place  immediately.  Thus  the  heat  release  rate  is  much  faster  than 
in  normal  diesel  engines  and  high  compression  ratios  are  not  required. 

Due  primarily  to  the  need  to  limit  the  scope  of  the  present  study 
coupled  with  the  fact  that  only  one  experimental  analysis  of  hypergolic 
combustion  has  been  conducted  to  date  [1],  the  thrust  of  the  work 
reported  herein  is  the  development  and  validation  of  a  working  computer 
code  for  use  in  future  engine  development.  An  analysis  of  these  engine 
concepts  using  the  code  developed  will  be  the  subject  of  follow  on 
research, 

1,2  Previous  Investigations  Leading  Up  to  the  Present  Analysis. 

A  review  of  the  past  decade  and  a  half  of  literature  pertaining  to 
the  numerical  modeling  of  combustion  in  IC  engines  and  on  related  topics 
reveals  that  there  are  a  maze  of  approaches  being  pursued  that  have 
applicability  to  this  problem.  Included  are  five  major  areas: 

1,  The  degree  of  geometric  dimension  used  in  the  analysis;  Zero 
dimensional  models  average  the  spacially  varying  quantities  and 
use  only  time  as  the  independent  variable.  One  dimensional 
models  use  time  along  with  only  one  special  coordinate  and 
involve  averaging  along  the  other  coordinates.  The  multi¬ 
dimensional  approach,  also  referred  to  as  full  simulation 


approach,  involves  the  simulation  of  all  processes  which  occur 
in  the  engine  including  fluid  motion. 

2.  The  engine  feature  simulations  considered;  This  includes  a 
variety  of  items  such  as  valve  models,  fuel  injection  models, 
spark  models,  flame  propagation  models,  wall  quenching,  piston 
motion,  etc. 

3.  The  fluid  dynamics  assumptions  made;  This  includes  such  items 
as  whether  the  flow  is  considered  viscid  or  inviscid,  laminar 
turbulent,  having  uniform  gas  properties  such  as  specific  heat 
and  diffusivity  or  properties  that  vary  spacially,  etc. 

4.  The  type  of  combustion  model  employed;  Turbulent  mixing  models 
and  chemical  kinetics  models  are  in  use.  The  chemical  kinetics 
models  vary  from  use  of  single  equation  models  to  use  of 
detailed  chemical  kinetics  which  considers  the  myriad  of  branch 
reactions  involved  in  hydrocarbon  combustion, 

5.  The  numerical  method  employed  to  solve  the  resultant  set  of 
model  equations;  The  approaches  in  use  vary  greatly. 

A  complete  accounting  of  the  references  pertaining  to  all  the  per¬ 
mutations  of  these  five  major  areas  would  be  excessive.  Instead,  a 
literature  route  which  briefly  describes  the  wake  of  significant  work 
behind  the  present  study  of  hypergolic  combustion  in  a  heat  barrier 
piston  engine  will  be  given. 

Work  begun  by  Blaser  and  added  to  by  Pouring,  Fail  la,  Rankin, 
et  al,  lead  to  the  development  of  an  innovative  internal  combustion 
engine  referred  to  as  the  "Heat  Balanced"  engine,  which  was  found  to 
give  efficiency  increases  of  over  40  percent  compared  to  conventional  IC 
engines  under  certain  operating  conditions.  The  major  feature  of  this 


design  is  the  carving  out  of  a  volume  from  the  piston  just  below  its  top 
and  along  its  outside  radius  creating,  in  effect,  a  second  chamber.  It 
was  found  that  with  this  design  a  portion  of  the  heat  which  normally 
would  have  been  lost  through  the  bottom  of  the  piston  was  instead  re¬ 
added  to  the  air  mixture  during  the  intake  and  compression  strokes  of 
the  combustion  cycle. 

The  discovery  of  this  phenomenon,  called  internal  regeneration,  led 
to  a  series  of  follow  on  studies.  Work  by  Adams^’^  to  study  heat 
transfer  in  the  "cap"  portion  of  this  piston  by  modeling  the  cap  as  a 
radial  fin  gave  insight  into  this  phenomenon.  A  theoretical  study  was 
conducted  by  Keating,  et  al,®’^  on  the  potential  effects  of  several 
assumed  mechanisms  on  engine  performance  using  the  air  standard  cycle  as 
the  basis  for  the  description  of  engine  stroke  processes.  The  model 
described  herein  is  the  result  of  a  parallel  effort  which  was  begun  at 
about  the  same  time.  It  was  developed  for  the  ultimate  purpose  ' 
studying  the  mechanism  of  internal  regeneration  and  represents  a  solu¬ 
tion  to  the  conjugate  heat  transfer  problem  between  the  piston  and  the 
engine  gases.  It  uses  a  full  simulation  approach  to  model  the  engine 
stroke  processes. 

In  the  search  for  an  adequate  engine  model  for  this  study,  it  was 
determined  advantageous  to  configure  the  engine  as  shown  in  figure  1-1, 
with  a  heat  barrier  placed  in  the  piston  a  distance  below  its  surface 
making  it  a  heat  barrier  piston  engine.  The  heat  barrier  piston  engine, 
a  type  of  adiabatic  piston  engine,  also  stores  and  re-adds  thermal 
energy  to  the  engine  cycle  and  is  considered  optimum  for  the  study  of 
this  concept  because  of  the  simplified  geometry  involved.  For  reasons 
explained  in  the  next  paragraph,  it  was  also  decided  best  to  configure 
the  engine  for  hypergolic  combustion.  There  has  been  much  recent 
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investigation  [10  -  18]  of  engine  modifications  similar  to  those  con¬ 
sidered  in  this  pi'esent  work.  The  thrust  of  the  research  dealing  with 
low  heat  rejection  (adiabatic)  engines  has  been  experimental  only 
[10  -  16]. 

A  key  to  the  study  of  cycle  processes  is  the  proper  modeling  of  spa 
cial  heat  release  throughout  the  combustion  stroke.  The  problems  with 
adequately  modeling  combustion  in  IC  engines  are  immense.  The  desire  to 
create  a  model  which  could  be  confirmed  by  future  experimental  work, 
coupled  with  the  desire  to  overcome  the  combustion  problems  associated 
with  the  overall  higher  operating  temperatures  of  the  heat  barrier 
engine  led  to  the  adoption  of  a  potentially  significant  combustion 
technique  envisioned  by  Hoppie^^  and  demonstrated  to  be  feasible  by 
Scharnweber  ’  .  The  adoption  of  this  technique  called  hypergolic  com¬ 

bustion  made  possible  the  use  of  the  fast  chemistry  model  described  by 
Lockwood  in  [20]  with  the  expectation  of  good  comparative  results  with 
future  experimental  work. 

The  advent  of  the  full  simulation  approach  is  fairly  recent. 
Previous  workers  modeled  engine  processes  using  zero  dimensional  [21,22] 
and  one  dimensional  [23  -  30]  approaches.  Though  much  work  along  these 
lines  continues  [31  -  35],  the  use  of  the  multi-dimensional  approach  is 
increasing  in  combustion  engine  analysis.  Significant  pioneering  work 
in  the  development  of  this  approach  was  conducted  by  groups  at  the 
University  of  Maryland  [36  -  40],  Imperial  College  [41  -  49],  Princeton 
[50  -  54],  Los  Alamos  [55  -  57],  Lawrence  Livermore  [50,52,58],  Science 
Application  Inc.  [59  -  61],  etc.  The  group  headed  by  Anderson  at  the 
University  of  Maryland  was  the  first  to  simulate  all  4-strokes  in  a 


firing  spark  ignition  engine  using  an  explicit  finite-difference  proce¬ 
dure  by  MacCormack^'^.  They  also  developed  the  first  3-dimensional 
engine  simulation.  Gosman,  et  al ,  at  Imperial  College  developed  a  code 
called  RPM  based  on  an  implicit,  iterative,  finite-difference  scheme 
and  has  used  the  code  to  perform  calculations  for  a  number  of  engine 
configurations.  The  work  by  Bracco,  Gupta,  et  al ,  at  Princeton  as  well 
as  Butler,  et  al ,  at  Los  Alamos  was  based  on  the  RICE  (Reactive  Implicit 
Continuous-Fluid  Eulerian)  technique  of  Rivard,  et  al,^^.  Work  by  Boni, 
et  al ,  of  Science  Application  used  a  semi-implicit  Arbritrary  Lagrangian 
Eulerian  (ALE)  method  of  Hirt,  Amsden,  et  al,^"^  in  calculations  for  com¬ 
bustion  in  a  strati fied-charge  engine. 

More  recently,  Haselman  of  Lawrence  Livermore  has  produced  a 
method  called  "TOC"  for  use  in  the  study  of  reacting  flow  using  an 
explicit  numerical  scheme.  Other  workers  of  note  in  the  field  include 
Ramos  and  Sirignano^^  '  who,  using  a  computational  scheme 
vary  similar  to  Gosman's,  carried  out  among  the  other  things  a  detailed 
analysis  of  flow  in  the  valve  region  of  a  motored  engine.  Bernard^^ 
applied  his  turbulence  closure  method  given  in  reference  [72]  to  examine 
non-reacting  turbulent  flow  in  the  compression  stroke  of  an  IC  engine. 

At  General  Motors,  Diwakar^^’^**  has  used  the  "CONCHAS  SPRAY"  method  of 
Butler,  et  al,^^  to  study  the  Direct- Injector  Stratif ied-Charge  (DISC) 
engine. 

Of  those  works  among  the  above  in  which  reacting  flows  was  con¬ 
sidered,  most  modeled  the  combustion  process  using  an  Arrhenius  kinetics 
formulation,  applicable  to  the  laminar  combustion  process.  Gosman, 
et  al,****"'*®  has  pursued  an  alternative  approach  explained  in 
reference  [20]  by  assuming  that  in  most  instances  turbulence  mixing,  and 
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not  chemical  kinetics,  dominates  the  local  combustion  rate.  In  the  case 
of  hypergolic  combustion,  this  assumption  would  be  even  more  valid  and 
is  thus  the  approach  adopted  in  the  present  analysis. 

The  emphasis  of  all  of  the  full  simulation  models  mentioned  above 
was  the  flow  field  of  the  engine  only.  No  consideration  was  given  to 
heat  conduction  inside  the  cylinder  walls  and  piston  of  the  engine. 

Thus,  prior  to  the  beginning  of  the  work  reported  on  herein  the  need  for 
modeling  this  conjugate  heat  transfer  problem  using  a  full  simulation 
approach  existed  in  order  to  study  the  flow  of  heat  in  low  heat  rejec¬ 
tion  engines. 

1.3  Scope  and  Contributions  of  the  Present  Investigation 

As  mentioned  in  section  1.1  the  objective  of  this  investigaion  is 
the  development  of  a  comprehensive  computer  code  for  use  in  the  study  of 
heat  transfer  in  an  IC  engine  configured  with  a  heat  barrier  inside  the 
piston  and  powered  by  fuel  which  is  sprayed  into  the  engine  at  a  tem¬ 
perature  above  the  fuel's  activation  temperature.  The  approach  taken  to 
the  work  has  been  to  adapt,  where  possible,  current  developments  in 
numerical  heat  transfer  [76],  turbulent  flow  modeling  [41,66,77,78], 
combu.Lion  modeling  for  the  case  of  fast  chemistry  [20,48,79],  turbulent 
boundary  layer  theory  [80],  etc,  in  the  formulation  of  the  computational 
model  in  order  to  develop  as  accurate  a  code  as  possible.  Compromises 
had  to  be  made  in  light  of  the  need  for  computational  efficiency.  Where 
no  available  developments  were  suitable  for  the  code,  original  modeling 
was  conducted.  Also  included  in  the  work  was  the  development  of  a  graphics 
package  for  use  in  displaying  and  analyzing  data  generated  by  the  code. 

The  following  major  assumptions  were  made  in  the  development  of  the 
model  (See  figure  1-1) : 
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gure  1-1.  Heat  Barrier  Piston  Engine 
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1)  Axisymmetric  fluid  flow. 

2)  Effects  of  swirl  neglected. 

3)  Effects  of  radiation  neglected  (not  a  good  assumption  if  the 
barrier  in  the  piston  strongly  retards  heat  flow  through  the 
bottom  of  the  piston).  The  barrier  used  in  the  present 
analysis  is  weak  and  thus  makes  this  assumption  reasonable. 

4)  Use  of  one  valve  placed  along  the  symmetry  axis  for  both  intake 
and  exhaust. 

5)  Use  of  law  of  the  wall  for  specification  of  momentum  and  energy 
equation  boundary  conditions.  (This  specification  only  really 
valid  for  steady,  1-D  flow  with  a  zero  pressure  gradient  in  the 
boundary  layer.) 

6)  Effects  of  wall  flame  quenching  neglected.  Higher  wall  tem¬ 
peratures  due  to  the  heat-barrier  make  this  a  reasonable 
assumpt ion. 

7)  The  time  scale  of  the  chemical  kinetics  is  considered  infinitely 
greater  than  the  time  scale  of  turbulent  mixing  and  thus  makes 

a  fast  chemistry  model  usable  throughout.  Note  that  because 
of  the  hypergolic  nature  of  the  combustion  this  is  believed  to 
be  what  is  actually  happening. 

8)  Specification  of  constant  temperatures  at  a  given  distance 

within  the  piston  and  cylinder  linings.  This  is  based  on 
experimental  findings  of  Eaton  Corporation  [81].  See  chapter  2 
for  the  detai 1 s , 

9)  Use  of  Perfect  Gas  Law. 

7  7 

10)  Use  of  k-e  turbulence  model  due  to  Jones  and  Launder  .  (This 
model  assumes  turbulent  kinetic  energy  is  isotropic  which  is 
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not  the  case  during  the  rapid  expansion  and  compression 
processes  taking  place  in  IC  engines.) 

The  above  assumptions  served  to  focus  the  development  of  the  code  toward 
a  very  specific  problem.  However,  the  necessary  considerations  taken 
into  account  in  the  development  of  this  code  made  its  scope  very  broad. 
For  this  reason,  to  limit  the  extent  of  the  present  work  to  within  reason 
an  indepth  analysis  of  this  special  engine  was  not  included.  It  will  be 
the  object  of  future  research. 

The  resultant  code  was  verified  for  accuracy  by  two  means.  The 
first  of  these,  explained  in  Chapter  7,  is  an  extension  of  an  error  ana¬ 
lysis  technique  developed  by  Shih*^^.  Details  of  the  results  of  this 
comparison  are  contained  in  Chapter  8.  The  second  verification  was  the 
comparison  of  the  results  of  a  computation  for  a  given  engine  con¬ 
figuration  with  those  of  another  researcher  in  the  field  whose  results 
for  this  same  computation  gave  good  correlation  with  experimental 
results  [69].  Details  of  this  comparison  are  also  given  in  Chapter  8. 

The  resultant  model  described  on  the  pages  that  follow  contains  a 
number  of  innovations  which  may  serve  to  advance  the  field  of  IC  engine 
research.  A  summary  of  these  contributions  listed  below  serves  not  only 
to  highlight  the  results  of  the  work,  but  also  to  add  clarification  to 
the  exact  scope  of  the  research  effort: 

1.  The  development  of  a  method  of  simultaneously  solving  heat  flow 
due  to  conduction  within  the  metal  of  the  cylinder  and  piston 
linings  and  heat  flow  within  the  engine's  gas  mixture: 

Previous  to  the  beginning  of  this  work,  researchers  in  the 
field  using  a  full  simulation  approach  specified  a  constant 


wall  tenperature  and  neglected  heat  conduction  within  the 


engine  wall  linings.  This  particular  development  makes 


possible  a  parametric  study  of  the  effect  of  varying  the 


strength  of  the  barrier  in  the  piston.  Also,  it  makes  possible 


the  study  of  the  effects  of  placing  barriers  at  other  locations 


in  the  engine.  In  order  to  complete  this  development  it  was 


necessary  to  solve  the  problem  of  how  to  transform  the  non-flow 


coordinate  system  for  the  energy  equation  to  facilitate  the 


matching  of  boundary  condition  nodal  points  with  the  trans¬ 


formed  flow  region  equations.  This  matching  was  required 


because  the  nodal  points  used  in  the  solution  of  the  governing 


flow  field  equations  are  moving  due  to  piston  motion. 


2.  The  development  of  a  governing  equation  transformation  tech¬ 


nique  which  enables  the  distances  of  the  near  wall  boundary 


nodal  points  to  the  walls  of  the  piston  and  cylinder  surfaces 


to  remain  unchanged  throughout  the  engine  cycle:  This  involved 


the  breaking  of  the  flow  fluid  into  three  solution  regions. 


3.  The  development  of  a  detailed  working  valve  model  which  enables 


the  analysis  of  heat  transfer  within  as  well  as  to  and  from  the 


valve  when  it  is  both  seated,  and  while  it  is  moving  through 


the  flow  field  region:  The  only  serious  valve  modeling  done 


prior  to  the  beginning  this  work  modeled  the  valve  as  infinitely 


thin  for  the  purpose  of  heat  transfer  and  neglected  energy 


equation  boundary  conditions  in  the  vicinity  of  the  open  valve 


[69]. 


4.  The  development  of  a  model  for  the  injection  process  of  fuel 


under  hypergolic  conditions:  Because  the  gaseous  fuel  at  the 
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injection  orifice  openings  was  determined  to  be  entering  the  com¬ 
bustion  chamber  at  Mach  1;  this  was  by  far  the  most  time  con¬ 
suming  and  challenging  of  the  obstacles  to  overcome  in  the 
solution  of  the  problem.  Involved  in  the  modeling  process  was 
the  development  of  a  special  discretization  scheme  for  the  nodal 
points  adjacent  to  the  fuel  entry  point.  Also  involved  was  the 
detailed  analysis  and  modification  of  governing  equation  source 
terms  in  the  vicinity  of  the  injection  point. 

The  validation  of  the  numerical  accuracy  of  the  power  law  discre¬ 
tization  scheme  for  the  convection  diffusion  problem  developed 
by  Patankar  :  An  extension  of  the  computational  scheme 
error  analysis  technique  developed  by  Shih  was  made  to  the 
computer  code  solution  of  the  five  applicable  discretized 
governing  equations  for  the  non-reacting  case.  This  analysis 
resulted  not  only  in  the  verification  of  the  accuracy  of  the 
code  but  also  pointed  to  the  preferred  utility  of  Patankar's 
power  law  discretization  scheme  in  the  discretization  of  the  total 
flux  terms  of  the  property  transport  equations.  See  section 
4.1,2,  chapter  7  and  chapter  3  for  details. 

The  determination  of  a  baseline  combustion  engine  model  for  use 
in  the  future  study  of  low  heat  rejection  engines:  Not  only  can 
this  model  be  used  to  study  the  specific  engine  innovations 
outlined  earlier,  but  because  of  its  simplistic  geometry  it  may 
also  be  used  as  the  basis  for  comparison  of  works  dealing  with 
more  complex  geometries,  such  as  the  heat  balanced  engine  [2-5] 
or  various  other  types  of  low  heat  rejection  engines  [10-16]. 

Also  it  can  be  used  as  a  basis  for  comparison  of  the  heat  flow 
characteristics  of  low  heat  rejection  engines  using  other  com¬ 


bustion  schemes. 


Chapter  2 

GOVERNING  EQUATIONS 

Calculation  of  heat  flow  in  the  fluid  flow  field  of  an  internal 
combustion  engine  is  governed  by  a  coupled  set  of  partial  differential 
equations  that  account  for  mass  continuity,  mass  momentum,  energy  and 
pertinent  chemical  species  mass  fractions  for  compressible  turbulent 
reacting  gas  flow.  Accounting  for  turbulence  using  the  k  -  e  model 
requires  the  addition  of  two  more  equations.  In  the  non-flow  portion  of 
the  engine  (i.e.,  cylinder  lining,  piston  and  valve)  only  the  energy 
equation  is  required. 

Piston  motion  within  the  engine  neccessitates  a  special  transfor¬ 
mation  of  the  governing  equations  to  facilitate  the  convenient  use  of  a 
control  volume  approach  in  the  numeric  solution  of  the  governing 
equations.  The  axisymmetric  geometry  of  the  engine  and  tne  neglecting 
of  the  effects  of  swirl  made  possible  the  solution  of  the  problem  in 
cylindrical  coordinates  using  only  the  r  and  z  directions  and  two  com¬ 
ponents  of  velocity. 

The  aim  of  the  present  chapter  is  to  delineate  the  governing  dif¬ 
ferential  equations  along  with  their  boundary  conditions  used  in  the 
formulation  of  the  problem.  Not  included  are  the  specifics  of  the  com¬ 
bustion  model,  and  all  the  boundary  conditions  associated  with  the  valve 
and  fuel  injector.  Also,  not  included  is  the  procedure  for  determining 
pressure.  These  items  will  be  covered  later.  The  governing  equations 
are  presented  in  both  their  untransformed  and  transformed  forms  for 
completeness . 
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Radial  Momentum  Equation  [41,  67] 


k-e  Turbulence  Model  Equations  [41,67,77,78] 

For  the  turbulent  kinetic  energy  (k): 

9  3  1  3  3  -  Upf^f 

—  (pk)  +  ^  (puk)  +  -  —  (rpvk)  =  —  — —  — 

at  az  r  ar  az  ai^  az 


^  [p  i^eff  3k 
r  ar  pj^  ar 

+  Si,  ,  (2-4) 


where 


Pk  =  1.00  , 


Sk  =  -Co  pe  +  G  , 

r  o  ^r^u-iz.ravio.vo  l  au  aV' 

5  =  2  (y  *  1-1^  1;)^  t 


7»_u  (  uef  f  +  pk ) 


Co  =  1.00. 


For  the  turbulent  kinetic  energy  dissipation  rate  (e): 

a  a  la  8  .ueff  1  a  ^  ueff 

--  (pe)  +  —  (pue)  + - (rpve)  =  —  [— - ]  + - [r  ^ - ] 

at  az  r  ar  az  pf-  az  r  ar  p^  ar 


where 


(2-5) 


Sg  =  —  (Cl  G  -  C2  p  e)  +  pe7*_u  , 
k 


(C2  -  Cl)  C.^2 


=  1.2  , 


Cl  =  1.44  , 
C2  =  1.92 


K  =  .4187. 


Energy  Equation  [66,67,77] 


,|(»T)  ♦  ,|(»l.T)  *i^(rpvT) 


3  r 

IT  1 


3Ti  ,  1  9  r  3T 

Wf  Jli  ^71F  —r 

+  St  ,  (2-6) 


where 


and 


V]  Mt 

Yeff  =  -  +  —  , 

Prt 

Pri  =  !•  , 

Prt  =  1*  . 

,  1  r^P  ap  ap  ^  *  ..•  *  ..'i 

Cp  at  az  3r 

•  ^ 

Qj-"  =  heat  addition  due  to  combustion  defined  in  Appendix  A  and 
modeled  in  Chapter  3, 

qr"  =  heat  addition  due  to  radiation  (neglected  in  present  analysis) 


Cp  =  specific  heat  (defined  in  Chapter  3). 


The  basis  for  formulation  of  equation  (2-6)  is  given  in  Appendix  A. 


Species  mass  fraction  and  mixture  fraction  [20,44,47,48] 
For  the  mass  fraction  of  i  (mi): 

1  a 

at  az  r  ar 


a  a  la  a  ueff  am^,  1  3  r  wpff  Sfiii  i 

—  (pmi)  +  —  (pumi)  + - (rpvmi)  =  -  [-- - + - [r  - - ] 

az  Sc  az  r  ar  ^  Sc  ar 


Smi  .  (2-7) 


where 


Sc  =  1.  (Schmidt  number)  , 


=  rate  of  species  generation  per  unit  volume  (modeled  in  Chapter  3). 


For  reasons  discussed  in  Chapter  3,  only  an  equation  for  the  mass  fraction 
of  the  fuel  (mf)  in  conjunction  with  an  equation  for  the  determination  of 
the  mixture  fraction  (f)  are  required  to  account  for  all  of  the  species 
concentrations  of  concern  in  this  problem. 

The  governing  equation  for  the  mixture  fraction  is: 


3 

at 


(pf)  +  —  (puf)  +  (rpvf) 

3z  r  3r 


-J.  [ilili  ll]  +  i  [r  il] 

32  Sc  3z  r  3r  Sc  3r 


(2-8) 


Thus,  for  this  governing  equation,  the  source  term  is: 

Sf  =  0  . 

Note  that  f  is  defined  and  explained  in  Chapter  3. 

General  Transport  Equation  [41,43,76] 

Equations  (2-2)  through  (2-8)  all  have  the  same  general  form  and  can 
thus  be  represented  compactly  by  a  more  general  transport  equation  in 
terms  of  the  primitive  variable  ji  as: 


3  3  13 

_  (p.p)  +  _  (pu4,)  +  -  —  (rpv(|)) 

3t  32  r  3r 


3  341 

3Z  ^  3Z^ 


13,  3$, 

^  — ]  +  S*  (2-9) 


r  3r 


3r  ■ 


where 


and 


4)  =  u,  V,  k,  e,  T,  mf,  or  f  , 

S,j,  =  corresponding  source  terms 

Fiji  =  corresponding  diffusion  coefficients. 


Adequate  in  internal  combustion  analysis  as  an  equation  of  state  is  the 


general i zated  state  equation  for  a  multi -component  ideal  gas  [83]: 


T  I  mi  R, 

i=l 

where 


(2-10 


Ri  =  — 


R  N-m 


Wi  Kg-K 


R  =  Universal  Gas  constant 


Wi  =  molecular  weight  of  species  i  . 
2.2  Transformation  of  flow  field  equations.  [67] 


Piston  motion  necessitates  the  transformation  of  the  z  coordinate  of 
the  computational  mesh  into  a  coordinate  system  which  does  not  change  with 
time.  In  order  to  ensure  optimum  treatment  of  the  boundary  conditions  at 
the  piston  and  cylinder  head  surfaces  however,  it  was  determined  pre¬ 
ferable  not  to  have  the  adjacent  nodal  points  change  in  their  relative 
distances  to  these  surfaces.  Further,  to  accommodate  valve  motion,  it  was 
determined  necessary  to  create  a  region  in  which  the  actual  mesh  point 
spacing  is  constant  in  time  in  the  z  direction.  These  three  problems 
generated  the  need  to  solve  the  governing  equations  in  three  separate 
solution  regions,  one  in  which  the  equations  are  not  transformed,  one  in 
which  the  spacing  of  the  computational  mesh  changes  with  time  in  accor¬ 
dance  with  piston  motion,  and  one  in  which  the  actual  spacing  between  grid 
points  in  the  z  direction  is  fixed  but  the  actual  positions  of  these  nodal 
points  are  changing  in  the  z  direction  with  piston  motion.  Further 


discussion  of  this  three  region  approach  is  contained  in  Chapter  5  along 
with  an  indepth  description  of  the  regions.  The  purpose  of  this  section 
is  to  describe  the  mathematical  formulation  of  the  transformation  only. 


Namely  - 


r,  t)  >  4.(x,  r,  f 


within  the  applicable  domain  of  transformation. 

Between  the  lines  x  =  0  and  x  =  1  on  the  geometric  representation  of 
our  engine  flow  field  solution  region  depicted  in  figure  2-1,  one  can 


see  that: 


z-distl 


where 


6p  =  Zp  -  distl  -  dist2 


(2-11) 


(2-12) 


Zp  i  distance  of  piston  face  from  the  z  origin.  Zp  is  a  function  of 
time  and  is  described  analytically  by  equation  (B-1)  in  Appendix 


Note  that 


3z  6p  3x 


3t  3t'  5r 


(2-13) 


(2-14) 


where 


d6p  dzp 

Up  =  - =  -  and  is  the  piston 

dt  dt 

velocity. 


Up  is  expressed  as  a  function  of  time  by  equation  (B-2)  in  Appendix  B 


n 


Letting 


TT  =  u  -  X  Up  (2-15) 

and  substituting  equations  (2-13)  through  (2-15)  into  equations  (2-1)  and 
(2-9)  results  in  the  following: 

Transformed  Continuity  Equation 

—  —  (pfip)  + - (pu)  + - (rpv)  =  0  .  (2-16) 

6p  3t  6p  3x  r  3r 


Transformed  General  Transport  Equation 


- (p5o<<))  + - (pu<(i)  +  —  (rpv({i) 

6p  6p  3x  r  Sr 


1  3  j. 

6p  3x  6p  3x 


—  (rrp 
3r 


3it) 

3r 


+  ,  (2-17) 


where 


is  the  transformed  source  term  for  a  given  i}>, 
2 • 3  Flow  Field  Equations  in  Transformed  Coordinates. 


The  application  of  the  coordinate  transformation  described  in 
section  2.2  to  the  specific  governing  flow  field  equations  of  our  problem 
was  simplified  by  the  fact  that  use  of  equation  (2-17)  made  necessary  only 
the  transformation  of  the  source  terms.  The  results  of  this  transformation 
are  given  in  Table  2-1.  Thus  Table  2-1  used  in  conjunction  with  equation 
(2-17)  gives  a  compact  summary  of  the  transformation  of  equations  (2-1) 
through  (2-8). 

2.4  Non-Flow  Energy  Equation 

The  analysis  of  heat  transfer  in  the  non-flow  portion  of  our  engine, 
consisting  of  the  piston,  valve  and  cylinder  lining,  involves  the  solution 
of  nnlv  one  eauation.  the  unsteady  heat  conduction  eauation: 


TABLE  2-1  Transformed  Governing  Equations  for  the  Flow  Field* 


wnere 


81 

r  K  — 
8r 


(2-18) 


8T  8  ,  3T  ,  18 

pC  —  =  —  K  —  + - 

8t  8z  ^  8z  r  8r 


C  =  specific  heat  (specified  as  465  J/Kg-K) 
and 

K  =  thermal  conductivity  (specified  as  50  J/Kg-m-s). 
The  density  p  is  specified  as  7800  Kg/m^.  Note  that  because  the  sides 
of  the  cylinder  do  not  move  with  piston  motion,  the  transformation 
employed  in  sections  2.2  and  2.3  does  not  apply  here.  Instead,  a  special 
technique  is  used  to  match  boundary  conditions  between  the  flow  and  non¬ 
flow  regions,  and  equation  (2-18)  is  solved  in  the  z  coordinate  system. 

For  details  of  this  special  procedure  see  section  2.5.4  and  Appendix  C. 

2 . 5  Boundary  Conditions 

The  general  expressions  for  the  flux  of  a  primitive  variable  (4))  at 
a  wall  boundary  is  given  by; 

Jw  (^)  =  -  ^  ^  ’w  >  (2-19) 

where 

41  =  can  be  either  u,  v,  T,  mf  or  f  in  our  case, 

Jw(  i{i)  =  f  lux  of  (()  at  wal  1 , 
r<)>w  ■  diffusion  coefficient  of  at  wall, 
and 

y  =  the  direction  perpendicular  to  the  wall  and  may  be  either  z  or  r 
in  our  solution  scheme. 

For  the  momentum  equations  (i.e.,  for  ()>  =  u  or  v),  r,j|  =  ugff  ^nd 
dw  "  ■'^w*  tde  shear  stress  of  the  wall.  In  the  case  of  the  energy  equation 

•  II 

(41  =  T).  =  ygff  and  =  q^  ,  the  rate  of  heat  flux  per  unit  area 

perpendicular  to  the  wall.  For  the  species  mass  fraction  equations 


.  .  r.-.'.i 


i 

.1 


V 

_  . 


-■1 


■i 


-  •  •  ( 
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In  the  discretization  of  the  general  transport  equations  to  be  carried  out 
in  section  4.1.2  using  a  finite  difference  approach,  the  flux  of  a  primitive 
variable  across  a  control  volume  boundary  will  include  the  sum  of  both  the  con¬ 
vection  and  diffusion  of  that  variable.  However,  at  the  wall  surfaces  within 
our  engine,  there  is  obviously  no  convection  and  the  finite  difference  approxi¬ 
mation  of  the  wall  flux  is  given  by:  [84] 


'^w  ( I* )  "  ■  I'd 


( tq  -  4'w) 


iSyw+g 


(2-20) 


where  (see  figure  2-2) 


(Jig  =  value  of  (j  at  the  flow  field  nodal  point  immediately  adjacent  to 
the  wal 1 , 

=  value  of  j)  at  the  wall, 

'5yw*g  =  distance  between  the  wall  and  the  adjacent  flow  field  nodal 
point 


=  value  of  diffusion  coefficient  at  the  wall  obtained  by  assuming 
that  the  law  of  the  wall  expres->ions  for  steady,  zero  pressure 
gradient  and  incompressible  flow  are  applicable, 
figure  2-3  depicts  the  law  of  the  wall  profile  for  velocity  used  in  this 
analysis.  Note  that  we  are  using  a  two  region  approach  by  defining  a 


r.  r;: 
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Figure  2-2.  Nornericlalure  for  Geometry  ol 
Near  Wall  Boundary  Coridilibris 


point  y  =  11.63  and  assuming  the  flow  to  be  purely  turbulent  above  this 
value  and  purely  viscous  below  this  point.  The  profile  for  temperature 
is  depicted  similarly. 

2.5.1  Velocity  Boundary  Conditions  [78j 

The  value  for  is  found  by  assuming  that  the  law  of  the  wall  rela¬ 
tionships  apply  for  our  case.  Close  to  the  wall  it  is  assumed  that: 

T  =  Tw  . 

The  quantity  y"*"  used  in  the  law  of  the  wall  formulation  is  evaluated  in 
this  analysis  by  using  the  value  of  the  sheer  stress  in  the  inertial  sub¬ 
layer  (tj)  for  Tw  so  that: 


-5yw*q  A  /- 


^  p  fiyw^g  Pg 


(2-21) 


where 

kg  =  value  of  turbulent  kinetic  energy  at  the  nodal  point 
adjacent  to  the  wall  (see  figure  2-2). 

This  follows  from  the  fact  that  in  the  inertial  sublayer  (30  <  y'*'  <  400) 
the  flow  is  assumed  to  be  fully  turbulent  and  the  local  rate  of  turbulence 
production  is  balanced  by  the  rate  of  viscous  dissipation  so  that: 

Ti  =  p  k 

"Strictly  speaking"  however,  the  expression  for  y'*'  given  in  equation 
(2-21)  above  is  only  valid  in  the  inertial  sublayer.  [78] 

Near  the  wall,  a  1-0  Couett  flow  analysis  (appropriate  for  such  con¬ 
ditions)  tells  us  that 


T  =  (ui  +  ut)  — 


For  y"*"  <  11.63  it  is  assumed  that: 


ui  >>  wt 


and 


(2-22) 


where 


Ug  =  velocity  parallel  to  the  wall  at  the  adjacent  nodal  point  to 
the  wal 1 . 

For  y'*'  >  11.63  it  is  assumed  that 

ut  >>  Ml 


and 

C^V4  pg  kg^2  K  Ug 

£  pg  kg  ^2 

loge  [ - 

Ml 


(2-23) 


where 


Pg  =  density  at  nodal  point  adjacent  to  wall, 

E  =  9.793  (value  of  integration  constant  based  on  smooth  wall  con¬ 
ditions) 


K  =  .4187. 

To  simplify  the  nomenclature  scheme  used  in  the  precise  presentation 
of  the  boundary  conditions  for  not  only  the  U  and  V  momentum  equation  but 
also  the  equations  for  k,  e,  T,  mf,  f  and  P  the  geometry  and  subscripts 
shown  in  figure  2-2  apply  with  the  following  modifications: 

For  the  cylinder  head  which  represents  the  west  (w)  side  of  the  solu¬ 


tion  scheme: 


Xwg  “  '^^w+g  > 
pgw  “  Pg  > 

+ 

and 

=  Jw  ('<’)  • 

For  the  cyclinder  wall  (side  of  cylinder)  which  represents  the  north 
(n)  side  of  the  solution  scheme: 

♦n  ■  C>w  • 

‘('gn  ~  'i’g  » 

•'ng  ■  '^yw*g  «  . 

Pgn  "  Pg  t 
+ 

rn  = 
and 

Jwp  (4')  ■  Jw  (4)  • 

For  the  piston  surface  which  represent  the  east  (e)  side  of  the  solu 
tion  scheme: 

4e  "  4w  > 

4ge  “  ‘t'g  j 

Xeg  =  'Syw-*^  » 

Pge  "  Pg  » 

+ 

xe  =  y^ 
and 


^Wp  (4)  =  Jw  (4)  • 


For  the  symmetry  axis  which  presents  tne  south  (s)  boundary  of  the 
solution  scheme  =  41  at  the  border  and  the  subscript  s  is  used  in  partial 


derivatives  to  indicate  south. 

The  above  nomenclature  scheme  is  generalized  for  4  and  (4)),  where 
Jyy  represents  the  appropriate  wall  flux  term  corresponding  to  a  particular 
primitive  variable  4i.  For  other  surfaces  such  as  the  valve  and  fuel 
injector,  similar  patterns  of  nomenclature  are  used.  The  boundary 


For  Tp  >  11.63: 


pgp  kgn^2  K  Uq 


loge  [E  rn 


Piston  Face 


Ug  -  Up  . 


<eg  Pge  '^ge^^ 


For  Xp  <  11.63: 


2-29) 


(2-30) 


Wl  Vge 


(2-31) 


For  xg  >  11.63: 


Pge  kgg^2  K  Vgg 


logg  [E  xg  ] 


(2-32) 


3u 

=  0  . 


(2-33) 


Vc  =  0  . 


Value  Surfaces  and  Orifice  Openinc 


See  Chapter  5. 


Fuel  Nozzle 


See  Chapter  6. 


2.5.2  k-e  Boundary  Conditions  [78] 


(2-34) 


As  mentioned  in  2.5.1  the  approach  adopted  here  is  strictly  valid 
only  for  the  inertial  sublayer,  where  the  assumption  is  made  that  the  flow 


iw; 


is  completely  turbulent  such  that  the  local  rate  of  production  of  tur¬ 
bulence  is  balanced  by  the  viscous  dissipation  rate.  Thus  the  value  of  e 
at  the  near  wall  nodes  is  fixed  as  follows: 


C,3/4  1,^3/? 


Obviously,  at  the  wall  the  values  of  e  and  k  are: 
=  0 


kw  =  0 


(2-35) 


(2-36) 


(2-37) 


There  is  no  contributory  flux  of  turbulent  kinetic  energy  from  the  wall 
for  the  k  equation  near  the  wall  and  integration  of  the  source  term  (S^) 
in  the  discretization  of  this  equation  requires  special  consideration  at 
this  boundary.  Discretization  of  Sk  at  locations  other  than  at  the  near 
wall  boundary  conditions  is  explained  in  Chapter  4.  Note  that: 


/  Sk  d( V)  =  -  /  (Cd  PG  -  "S)  d( V) 

V  V 


The  integration  of  e  is  modeled  as  follows  in  order  to  extend  the  k- 
balance  formulation  of  equation  (2-35)  into  the  buffer  area  in  addition  to 
the  viscous  sublayer; 


Cd  Pg  Cy3/4  kg3/2  v 

/  Cope  d(V)  =  - - -  U+ 

V  fiyw+g 


(2-38) 


where 


for  y"*"  <  11.63  - 


:  ¥  ft 


and  for  y"*"  >  11.63  - 


U"  =  i  loge  [E  y+]  . 

The  integration  of  G  is  varied  from  the  discretization  conducted  in 
Section  4.1.3  by  noting  that  at  the  wall  the  fourth  term  in  the  expression 
for  G  given  in  Table  2-1  can  be  integrated  as: 


3U  1  3V  p  (Ug  -  Uvv) 

/  ueff  (—  *  T  ^ - T -  ^  • 

y  3r  5p  3x 


(2-39) 


Following  the  same  pattern  given  in  2.5.1  in  terms  of  nomenclature, 
U'*’  is  expressed  as  follows  for  the  various  internal  surfaces: 

+ 

Vw  =  U'*'  on  the  cylinder  head, 

+ 

Un  =  U'*’  on  the  cylinder  wall 
and 


Vg  =  U'*'  on  the  piston  surface. 
Cylinder  Head 
+ 

For  Xw  <  11.63: 

+  + 

Vw  ■  *w  • 

+ 

For  xw  >  11.63: 

+  1  + 

Vw  =  —  logg  [E  Xy^l  . 


C  3/4  ^^3/2 
''■gw 

^5'^  '  K  X 

^wg 

(2-41) 

o 

II 

5 

(2-42) 

kgv<  is  determined  by  using  the  following  expression  for  /  S|(  d(V)  at  that 

V 

nodal  point: 

_  ^0  Pgw  Cp3/4  kgy^3/2  +  ^ 

I  Su  d(V)  = - ^ -  Vw  V  +  /  G  d(V) 

V  y  Xwg  V 

,  1  3v  2  Vgw 

-  I  ueff  (6p  T7  )  d{V)  +  ~x“  ^ 

(2-43) 

where 

/  "S'  d(V)  specified  in  4.1,3. 

V 

Cylinder  Wal 1  (Side) 

+ 

For  rp  <  11-63: 

+  + 

Un  ■  '"n 

+ 

For  rp  >  11,63: 

+  1  + 
un  =  -  loge  [E  rp]  , 

K 

Cn  =  0 

(2-44) 

C  3/4  kpp3/2 

kgn  determined  from; 


/  V  ow 

V 


Cp  Pgn 

■^ng 


3/2  4- 

-  un  V  +  /  b  d(V) 

V 


r  3u  2  ^wn  ^gn 

■'  ^'eff  (—  )  d(V)  +  - V  . 

V  3r  Xpg 


Pi ston  Face 
+ 

For  Xe  <  11.63; 

+  + 

Vg  =  Xg  . 

+ 

For  Xg  >  11.63; 

■*■1  + 

Vg  =  logg  [E  Xg]  . 
K 


eg  =  0  . 
ege  =  0  • 
k  g  -  0  . 

kgg  determined  from; 


/  ^kgg  (j(V) 

V 


Cd  Pge  Cy3/4  kgg 


3/2  vt 


xeg 


+  /  G  d(V) 
V 


/  i^eff  (■  T"/ 


d(V) 


■^Wg  Vge 


I  hnc  d(V)  =  -(C[)P£  -  V  -  /  Ueff  ( — )  d(V) 


(2-54) 


Value  Surfaces  and  Orifice  Qpeninc 
See  Chapter  5. 


Fuel  Nozzle 


See  Chapter  6. 


2.5.3  Temperature  Boundary  Conditions 

The  boundary  conditions  for  the  energy  equation  are  specified  at  a 
fixed  distance  within  the  interior  metal  surfaces  from  the  borders  of  the 
flow  field  solution  region.  Work,  done  by  Eaton  Corporation  [81] 
indicates  that  for  an  IC  engine,  at  a  specific  distance  6  inside  the  metal 
of  the  cylinder  or  piston  from  the  combustion  chamber,  the  temperature 
will  not  vary  over  1%  during  the  engine  cycle.  This  distance  varies 
depending  on  the  thermal  conductivity  of  the  metal  and  the  RPM  of  the 
engine.  For  the  cylinder  wall  (side  of  cylinder)  it  also  varies  with  the 
distance  from  the  cylinder  head  (i.e.,  varies  in  the  z  direction).  The 
actual  thickness  and  values  of  the  corresponding  temperatures  can  be 
determined  experimentally.  Because  experimental  data  does  not  exist  for 
our  engine,  assumed  values  are  used.  See  figure  2-4  for  the  values  of 
the  temperatures  specified  for  this  analysis.  The  values  of  the  various 
penetration  6's  are  as  follows:  6^  =  .0004  m;  63  =  .0003  m;  and  63  =  .0005  m 
The  other  two  6's  shown  in  figure  2-4  are  used  for  specification  of  the 
limits  of  the  temperature  borders  and  are  as  follows:  64  =  .03735  m; 

65  *  .05065  m.  Note  that  these  values  will  also  change  for  the  heat 


barrier  piston  engine  with  the  degree  of  barrier  which  is  placed  within 
the  piston. 

The  temperature  boundary  conditions  at  the  perimeter  of  the  fluid 
flow  field  are  thus  not  specified.  The  heat  flux  is  determined  at  this 
border  using  the  same  procedure  as  was  used  in  the  solution  of  the 
momentum-transport  equations.  The  wall  fluid  boundary  layer  is  assumed  to 
have  a  constant  heat  flux  throughout  so  that: 

•  •  II 

d  ~  dw  1 

where  for  1-D  analysis  (Couette  flow) 

3T 

d  =  (ti  +  Yt)  Cp  —  . 

3y 

Note  that: 


wl 


and 

wt 

Tt  =  -  • 

Prt 

Using  Reynolds'  analogy  and  assuming  that  the  viscous  sublayer  and  thermal 
viscous  sublayer  are  the  same  thickness,  for  y"*"  <  11.63  it  is  assumed 
that  [66,78,80]: 


Y1  »  Yt 


and 


*  II 


Wl 


Pr] 


(Tq  -  Tw) 


Also,  for  y"*"  >  11.63 


(2-55) 


Yt  >>  Y1 


(2-56) 


*  ^9  "^Pg  ^9^^  ^"^9  ' 

q"  - - 

w  Prt 

— ^  loge  (E  y+)  +  Pr.  (P  4  — ^  4-) 

K  Pr^ 


where  [80] 


o*  4 

P  4  —  f  =  - 

a*  t  ^ 

’  Sin  (-) 

4 


(-)V2  (  ^ 

K  Oa  ^ 


-  1)(  AL  )V4  . 


A  =  Van  Driest  Constant 
=  26.0  for  smooth  wall 


a<{,  =  Pr,  Sc,  etc,,  depending  on  4> 


Thus  for  the  various  inner  metal  surfaces,  the  flow  field  and  non-flow 
energy  equation  heat  flux  boundary  conditions  can  be  specified; 
Cylinder  Head 


For  <  11.63: 


'^W  Xy^g  Pp^ 


(^qw  -  Ty,) 


(2-57) 


For  Xyy  >  11.63: 


q  Ww  ■  ■ 


Pgw  "^pgyy  *^gw^^  (^gw  -  ^w) 

Pr*.  +  Pri 

— ^  loge  (E  xw  )  +  Prt  (P  4  —  f  ) 

K  Prt 


(2-58) 


^linden  Wall  (Side) 


For  rn  <  11.63: 


Piston,  Cylinder  Interface 


Heat  flux  between  the  nodal  points  of  the  side  of  the  piston  (ps)  and  the 
side  of  the  cylinder  wall  (cw)  at  the  interface  shown  in  figure  2-5  is 
modeled  as: 


'l"p-K:  "  ^^p+c  (Tps  "  T^cw)  , 


:2-64) 


where 


*^p+c  "  combined  convective  and  radiative  heat  transfer 

coefficient  between  the  piston  side  and  cylinder  wall.  It 
must  be  determined  experimentally  and  is  a  function  of  RPM, 
pressure,  temperature,  etc.  For  our  analysis  a  fixed  value 
was  selected. 

Note  that  the  value  of  Tp^  varies  along  the  left  portion  of  the 
piston  cylinder  interface  and  is  a  fixed  value  (T[  -  See  figure  2-4) 
along  the  right  portion  of  the  interface.  Obviously  Tpn,  and  the  left  por¬ 
tion  of  Tps  vary  with  position  and  time  and  are  computed  iteratively  along 

« 

with  q"p>c  in  the  solution  of  the  non-flow  energy  equations.  As  can  be 
seen  in  figure  2-5,  below  the  piston,  a  portion  of  the  cylinder  lining 
considered  in  the  non  flow  solution  region  is  bordered  by  air.  At  these 
nodal  points,  the  boundary  conditions  are  specified  by  equation  (2-59). 
Valve,  Valve  Seat  Interface 

A  detailed  description  of  heat  transfer  within  the  valve  is  given  in 
Chapter  5  as  previously  indicated.  Heat  flux  at  the  interface  of  the 


valve  seat  (vs)  and  the  closed  valve  is  given  by 
''side'*’''^  ■  ^Vsi(je>vs  -  Tvs)  » 


(2-65) 


In  l «  2.^. 


L  i 


I  -i 


t-.  J 


'  o  ..tjl 
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-  area  included  In  nonflow  solution 


Figure  2-5.  Nonflow  Solution  Region. 


where 


h^sicie*''^  =  comDined  convective  and  radiative  heat  transfer  coef¬ 
ficient  between  the  valve  side  and  the  valve  seat. 

Comments  for  equation  (2-64)  for  hp^  apply  here  also. 

When  the  valve  is  open,  air  flow  through  the  valve  may  either  be  into  or 
out  of  the  engine.  When  the  flow  is  into  the  engine 

^  rvs»g  Pin  *^in^^ 

rvs  = - - -  . 

Pin 

where 

Cy  Pin  l^^in 

Pin  =  Pi  +  - 

^i  n 

and 

Pin»  i^in  ^in  -  3re  the  near  wall  inlet  values  of  p,  k,  and  e 
respectively  as  specified  in  Chapter  5. 

Heat  flux  from  the  adjacent  gas  to  the  valve  seat  is  as  follows: 

+ 

f^or  rvs  11-63: 

'’"gvs^vs  =  srv^^ 

where 

Tg^n  "  temperature  of  atmosphere. 

+ 

For  rvs  >  11.63: 


(2-67) 


Q  gvs+vs 


Pin  ^Pir\  ^^9vs  ‘ 

P  +  P 

—  loge  (Ervs)  ^  Pr.  (P  4  —  f  ) 
K  P,, 


If  the  flow  of  air  is  out  of  the  engine,  the  substitution  of  the  near 
wall  out  flow  values  of  p  ,  k,  y,  e,  T  are  made  for  the  corresponding  in 
flow  values  in  the  above  equations  (equation  (2-66)  -  (2-67))  to  obtain 
equivalent  expressions  for  q"g^^*v5  in  the  out  flow  case.  See  Chapter  5 
for  a  specification  of  these  values. 


2.5.4  Fuel  Species  and  Mixture  Fraction  Boundary  Conditions 


The  flow  can  be  assumed  to  be  non-reacting  at  and  very  near  the  walls 
due  to  the  cooling  effects  of  the  wall.  This  is  normally  true  for  the 
piston  surface  also,  assuming  the  barrier  used  in  the  piston  is  such  that 
the  surface  temperatures  remain  cool  enough.  Under  such  conditions,  the 
specification  of  the  flux  of  mf  and  f  can  be  specified  in  a  manner 
identical  to  that  used  for  the  energy  equation  [80].  However,  to  use 
these  relationships,  the  values  of  mf  and  f  must  be  known  at  the  wall 
throughout  the  engine  cycle.  A  simple  approach  is  to  assume  that  mf  and  f 
are  botn  zero  at  the  walls  and  to  proceed  exactly  as  for  the  case  of  the 
energy  equation.  In  this  approach  Jfuel  and  Jf  are  found  by  the  use  of 
identical  equations  to  equations  (2-57)  through  (2-63),  by  using  Sc  in 
place  of  Pr  and  mf  or  f  in  place  of  T  to  find  Jfgel  Jf  place  of  q"y^. 
Initial  analysis  utilized  this  formulation.  However,  neither  mf  or  f  are 
zero  at  all  wall  locations  throughout  the  entire  engine  cycle.  If  their 
actual  values  could  be  known,  a  solution  using  this  formulation  would  be 
optimum.  Later  in  the  analysis  it  was  determined  better  to  assume  that 

AA 


the  values  of  mf  and  f  are  the  same  at  the  near  wall  and  wall  nodal 
points.  This  assumption,  though  not  optimum,  was  found  to  lend  itself  to  a 
much  more  accurate  solution  and  is  used  in  the  analysis  presented  herein. 
Cylinder  Head 
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'  w  'gw 

Jfw  =  0  . 

^w  "  fgw  • 
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See  Chapter  5. 


Fuel  Nozzle 


See  Chapter  6. 
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Chapter  3 


COMBUSTION  MODEL 

In  the  engine  model  under  investigation,  the  rate  of  combustion  is 
controlled  by  turbulent  mixing  and  not  chemical  kinetics.  This  is  a 
consequence  of  the  fact  that  the  fuel  enters  the  engine  at  a  temperature 
above  the  fuel's  reaction  activation  temperature.  Under  such  conditions 
(hypergolic  conditions),  relative  to  the  rate  of  turbulent  mixing,  the 
rate  of  chemical  reaction  is  so  fast,  it  can  be  viewed  as  occurring 
instantaneously  once  a  radical  of  fuel  and  oxygen  are  mixed  [1,18,19,20]. 
As  a  consequence,  the  array  of  intermediate  chemical  reactions  which 
occur  in  the  combustion  process  can  be  ignored.  The  only  reaction 
equation  of  concern  is  the  overall  heat -releasing  reaction  given  by: 

1  Kg  fuel  +  s  Kg  02  ♦  (1+s)  Kg  product.  (3-1) 

Due  to  its  high  injection  temperature,  the  fuel  is  assumed  to  react 
irreversibly  with  the  oxygen. 

The  effect  of  the  formation  of  nitrogen  oxides  during  the  com¬ 
bustion  process  have  relatively  little  effect  on  the  global  heat  release 
process  within  the  engine  and  as  a  consequence  is  neglected  [47]. 

Also,  to  simplify  analysis,  the  mass  diffusion  coefficients  of  each  spe¬ 
cies  are  assumed  equal.  This  last  assumption  is  well  satisfied  for  the 
case  of  fully  turbulent  flow  [20]. 

As  a  consequence  of  the  above,  our  combustion  model,  in  addition  to 
the  local  rate  of  fuel  combustion  (-mf),  need  only  be  concerned  with  the 
computation  of  the  local  mass  fractions  of  the  fuel  (mf),  oxygen  (mo2), 
and  combined  combustion  products  (mpr).  With  the  rate  of  fuel  com- 
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bustion  specified,  the  rate  of  heat  release  per  unit  volume  q^"  can  oe 
determined  for  use  in  the  solution  of  the  energy  equation.  The  sections 
which  follow  delineate  this  model. 


3.1  Species  Mass  Fraction  and  Mixture  Fractions 

The  general  species  mass  fraction  equation  (given  in  its  untrans¬ 
formed  form  in  equation  (2-7)),  need  only  be  solved  for  one  species  by 

observing  first  of  all  that  a  passive  scalar  f,  defined  as 

"^02  mo2,a 

rnf  -  —  +  -T-  (3-2) 

f  ■  — rr^oxa— 

obeys  the  transformed  general  transport  equation  with  its  source  term 
(Sp)  equal  to  zero,  and  by  noting  that  given  f  and  one  species  mass  frac 
tion,  all  the  otner  mass  fractions  can  be  determined  algebraically.  In 
the  above  expression  mo2.a  is  mass  fraction  of  O2  in  air  under 
atmospheric  conditions  and  s  is  defined  by  equation  (3-1).  Once  f  is 
determined,  the  mass  fraction  of  oxygen  (mo2)  can  be  determined.  The 
mass  fraction  of  inert  gas  (mip)  is  obtained  from  the  fact  that  f  and 
min  ii3ve  the  simple  linear  relationship  given  by  [20,48]; 

'^in  “  '^in,a  (1“^)  (3-3) 

where  min^a  is  the  mass  fraction  of  inerts  in  air. 

The  mass  fraction  of  the  fuel  (mf)  is  obtained  by  solving  the 
transformed  version  of  equation  (2-7)  for  mf  (see  Table  2-1).  Needed  in 


this  calculation  is  a  model  for  mf,  the  rate  of  change  of  the  mass  frac¬ 
tion  of  fuel  per  unit  volume.  An  eddy-dissipation  combustion  model 
developed  by  Mangussen  and  Hjertager  for  a  single  step  irreversible 
global  reaction  is  given  by: 

•  £  ''*0?  ^  '’’pr 

mf  =  -A  (-)  min  (mf,  -  ,  - ),  (3-4) 

K  s  1  +  s 


where 

A  =  4  (obtained  exper imental ly) , 

B  =  .5  (obtained  experimentally), 
and 

min(,,)  =  minimum  numerical  value  of  the  expressions  contained  in 
parenthesis. 

This  is  considered  to  be  a  reasenable  model  for  hydrocarbon  fuels 
at  engine  conditions  and  its  utility  in  internal  combustion  engine 
research  is  demonstrated  in  reference  [48].  The  values  for  A  and  B  may 
require  some  adjustments  in  order  to  further  refine  this  expression 
(equation  (3-4))  for  engine  use.  Also,  this  model  is  not  applicable  in 
the  case  of  either  an  excessively  lean  or  rich  mixture. 

The  final  mass  fraction  needed  for  our  solution  is  that  of  the  com¬ 
bined  products  (mpf)  which  is  calculated  by  use  of  the  fact  that  the  sum 
of  the  mass  fractions  must  equal  one. 

nipr  =  1  -  mf  -  mo2  “min.  (3-5) 

3 . 2  Energy  Equation  Combustion  Source  Term  Model 

As  shown  in  Appendix  A,  the  untransformed  form  of  q(-  •  tne 

contribution  of  combustion  to  the  energy  equation's  source  term,  is 
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given  by  the  expression 


qc 


n  .  0 

■  "  p  I  in-j  hi 
i  =  l 


(3-6) 


when  the  Prandtl  and  Schmidt  (and  thus  Lewis)  numbers  are  equal  to  one. 
In  our  case,  "i"  in  the  above  expression  can  represent  either  fuel,  oxy¬ 
gen  or  products. 

The  rate  of  production  of  fuel  per  unit  volume  (mf)  is  given  in  the 


preceding  section.  Equation  (3-1)  dicates  that 


mo2  =  s  mf 


and 


mpr  =  -  (1  +  s)  mf. 


(3-7) 


(3-8) 


0  0  0 

The  tabular  values  ho2.  hf  and  hp^  are  given  in  Table  3-1  for  the  use 

0 

of  C12H26  as  a  fuel.  Note  that  ho2  equals  zero.  Thus  our  expression 

*  II  ' 

for  q,.  reduces  to 

.0  0 

qc  =  -P  mf  (hp  -  (1  +  s)  hpr).  (3.9) 

3.3  Variation  of  the  Specific  Heat  (Cp)  Due  to  Combustion. 

In  the  formulation  of  an  energy  equation  based  on  temperature 
instead  of  enthalpy  (see  appendix  A),  the  specific  heat  of  the  gas  at 
constant  pressure  (Cp)  was  introduced.  Cp  varies  with  both  temperature 
and  species  concentration  in  accordance  with  the  relationship 


TABLE  3-1  Heats  of  Formation  [85] 


Species  (i 
Cl2  ^26 
02 

C02(g) 
ri20(g) 

Product 


Coversion 


Ci2  H26  +  I8V2  O2  ♦  12  CO2  +  13  H2O 


) 


0 

Heat  of  Formation  (hi) 


-  71,014  cal/  g  mole 
0 


-  94,051  cal/  g  mole 

-  57,798  cal/  g  mole 


Molecular  Weight  (Wi ) 
170.34 
31.9988 
44.0099 
18.016 


0  0 

12  hco;,  -*•  13  hHpQ 

25 


75,199.4 


cal 

g  mole 


30.49 


Factor; 


'i  (r-) 


Kg 


0  cal 

(hi  ( - --))  X  { 

g  mole 
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NUMERICAL  PROCEDURES 


The  numerical  method  employed  in  the  solution  of  the  governing 
equations  is  a  finite  difference  scheme  based  on  the  Teach  -  T  code  [78] 
with  major  modifications.  The  code  was  converted  from  a  self-selecting 
hybrid  central/upwind  spacial  differencing  scheme  to  a  power  law  spacial 
differencing  scheme  [76].  The  latter  scheme  uses  an  expression  which 
almost  exactly  approximates  the  solution  of  the  steady  one-dimensional 
convection  diffusion  problem  in  the  discretization  of  total  flux  across 
each  of  the  surfaces  of  the  control  volume.  The  code  was  modified  from 
one  which  handled  only  steady  non-reacting  flow  to  one  which  handles 
unsteady  flow  with  combustion.  The  source  terms  of  the  governing 
equations  used  in  the  codes  for  u,  v,  k,  e,  and  T  were  significantly 
modified  to  reflect  the  more  extensive  source  terms  listed  in  Table  2-1. 
Additionally,  the  code  was  only  usable  for  uniform  spacing  in  the  r  direc¬ 
tion,  used  coarse  averaging  techniques  for  viscosity  at  the  border  of  the 
staggered  velocity  control  volumes,  contained  errors  in  the  computation  of 
some  of  the  source  terms,  etc.  These  problems  were  resolved  through  code 
modifications. 

Details  of  the  discretization  and  solution  of  the  governing  equations 
are  given  in  the  succeeding  subsections  of  this  chapter,  including  the 
procedures  used  in  the  computation  of  pressure  throughout  the  flowfield. 

4.1  Discretization  of  Governing  Equations 

The  finite-difference  method  employed  in  this  analysis  uses  a 
staggered  grid.  In  this  formulation,  the  majority  of  the  dependent 
variables  are  calculated  at  the  center  point  (P)  of  the  control  volume. 

The  u  and  v  velocities  are  displaced  so  as  to  lie  ha  f  way  between  the 


pressures  that  drive  them  and  are  thus  on  the  faces  of  the  control  volumes 


used  in  the  calculations  of  the  other  primitive  variables.  See  figure 


4-1.  The  control  volimes  used  in  the  finite-difference  formulation  of  the 


momentum  equations  for  u  and  v  are  shown  in  figures  4-2  and  4-3  respec¬ 


tively. 


4.1.1.  Discretization  of  the  Transformed  Continuity  Equation 


The  transformed  continuity  equation  given  in  equation  (2-16)  is 


discretized  by  integrating  each  term  over  the  control  volume  shown  in 


figure  4-1.  Refer  to  this  figure  for  a  definition  of  the  geometrical 


terms  and  subscript  nomenclature  that  follow. 
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0  i  previous  time  step 
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.n  .n 
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and 

■  '^p  ‘^^ew  • 

Adding  the  above  terms  results  in  the  discretized  form  of  the  con¬ 
tinuity  equation: 


n  0  .n  .n  .n  .n 
Mp  -  Mp  +  m^  -  m^  +  mg  -  =  0  . 


(4-1) 


4.1.2  Discretization  of  the  Transformed  General  Transport  Equation 

Equation  (2-17)  is  integrated  over  a  general  control  volume  after  com¬ 
bining  several  of  its  terms.  If  41  =  k,  e,  T,  mf,  or  f,  the  geometry  and 
subscript  nomenclature  appropriate  to  figure  4-1  is  used.  If  4  =  u,  the 
geometry  and  subscript  nomenclature  for  figure  4-2a  is  used  instead.  If  4) 
=  V,  that  for  figure  4-2b  is  used. 

13  [(p  V)p.  -  (P  ^  V)?  ] 

J  J  —  —  (P'5p4>)  fipr  3r  3x  - - 

V4,  6p  3t  At 


n  n 

/  /  —  —  (pu4i - i  —  )  6nr  3r  3x 

Va  '^p  3x  6p  3x 


0  0 


ly/  S(|,  d(V((,)  -  Sp^  (|>p  +  Su^  • 

•  A 


where 

SP(|,,  Su^  -  are  specified  for  various  i|)  in  section  4.1.3 
Collecting  terms  yields  - 
n  n  0  0 

-  Js^  =  <t>p  +  Su^  , 

Note  that  if  equation  (4-1)  is  re-derived  for  a  generalized  control 

volume  in  a  manner  identical  to  the  formulation  of  the  above  equation,  and 

the  result  is  first  multiplied  by  (jip  and  then  subtracted  from  the  above 

w 

equation,  the  following  expression  results: 

0  n  0  .n  n  .n  n 

^P^,  ■  '<’P4,)  ^  'I'p^)  -  (Jw^  -  <j>p^)  + 

.n  n  .n  n 

("•n^  -  f>in^  4ip^  -  (Js^  -  ’’'s^)  ')>p^  )  =  Sp^  4.p^  +  Su^,  . 

The  above  result  is  simplified  by  using  the  power  law  spatial  dif- 

7  A 

ferencing  scheme  developed  by  Patanker  to  approximate  the  total  flux  terms 

•  n 

(0  -  m^  (|»p)^  jhg  total  flux  terms  become: 

.n  n 

-  "^e^  "kp^)  =  («j)p^  -  , 

.n  n 

(Jw^  -  ’^w^  <l'p^)  =  -  4ip^)  , 

.n  n 

('^n^  -  f^n^  tp<(,)  =  (4ip^  -  (|.n^) 

and 

.n  n 

(Js^  -  ms  <(»p^  )  =  as^  (<^5^  -  4>p^)  . 

Substitution  of  the  above  four  expressions  into  our  previous  equation 

results  in  the  final  version  of  the  discretized  general  transformed 

transport  equation  in  compact  form: 

no  n  n  o  o 

(3p^  Mp^  -  Sp^)  ♦p^  =  J  am^  4ni  +  Su^  +  Mp^  (4-2) 
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where  [76] 


^  j 

=  aE^  +  aw^  +  +  as^  ,  (4-2a) 

=  De^  A  (  I  PeJ  )  +  max  (-Ce^.  0)  .  (4-2b) 

^^4,  =  %  ^  U  %  I  )  0)  •  (^-2c) 

=  Dn^  A  (  I  Pn  J  )  +  max  (-Cp^.  0)  ,  (4-2d) 

=  Ds^  A  (  I  Ps  J  )  +  "lax  (Cs^,  0)  .  (4-2e) 


(4-2f) 


and 


A  (  I  P  I  )  =  max  (0,  (l-.l  I  P  I  )5) 


max  (,)  =  maximum  positive  value  of  quantities  contained  inside 
parenthesis 


For  clarity  it  should  be  noted  that  the  generalized  expression  for  the 
discretized  transport  equation  given  above  (equation  (4-2))  applies  for 
all  For  4i  other  than  u  or  v,  expressions  (4-2a)  through  (4-2i)  are 
derived  using  the  control  volume  shown  in  figure  4-1  and  the  nomenclature 

is  exactly  as  given  in  these  expressions  with  the  iji's  removed.  For  ^  =  u 

.n 

the  terms  a^^  (M  =  N,S,E,W),  Dhi^j,  c^y,  moy,  and  Pmy  (for  m  =  e,w,n,s)  are 
derived  instead  using  the  geometry  given  in  figure  4-2a.  Similarly,  for  41 


=  V,  the  terms  aiv^^,  Df^^,  rrip^  and  are  used  in  expressions  (4-2a) 
through  {4-2i)  and  are  defined  in  terms  of  the  geometry  given  in  figure 
4-2b. 

4.1.3  Discretization  of  Transformed  Source  Terms 

The  source  terms  of  the  transformed  governing  equations  (summarized 
in  Table  2-1)  are  discretized  by  integrating  each  over  its  appropriate 
control  volume. 

u-Momentum  Equation  Source  Term 

The  expression  for  the  transformed  u-velocity  transport  equation 
source  term  (Sj)  repeated  below  is  integrated  a  term  at  a  time  over  the 
control  volume  given  in  Figure  4-2a. 


.  J-iP +J:-i  [iielf  3ii]  +  i-i  [LJiaLL^]  .  JL  i  7.u  +  pk) 


6p  ax  6p  3x  6p  ax  r  3r 


where 


1  ap 

USl  i  /  /  (- -  }  5p  r  ar  ax 

«p  3X  ^ 

Vu 


1  (Pp  -  Pw) 

5p  6xpw 


'u  » 


(4-3) 


Vu  =  4p  rj  6rns  Sxewy  • 

US2  i  /  /  (^^  ])  d(v„) 

vp  9x  6p  3x 
Vu 


ueff 


weff  3u) 


where 


VS5 


Vy  =  6p  dxgvj  . 

VS2  H  /  J  (— —  [ueff  iii])  6p  r  3r  3x 

3x  3r 


r 

r  3u. 

li'eff  ~Je  ■ 

l-eff  ^1' 

Uy  '  Il5y 


VS3  =  I  J  (i  —  [r  Ueff  — ])  6p  r  3r  3x 
r  3r  3r  ^ 


3v  3v 

(r  I'eff  "^'ny  -  [>'  Ueff  "^JSy)  <5p  6xew 


3r^=‘v' 


VS4  =-  /  /  (-2  Ueff  — )  d  (Vy) 
Vy 


=  -2  ueffp^  Vy 


VS4'  = 


VS4 


,  r  ,  2  3  ,  ^  1  3u  V  3v , 

VS5  .  n  (-  3  -  (».ff  t-;^  ^  (“v) 


1  9li  1  3u 

(ueff  (t*”  ■*■  pk))p  -  (ueff  pk))s 


6p  3x 


6|^  3x 


6rps 


(4-10) 


(4-11) 


(4-12) 


(4-13) 


63 


kS3  =  j  J  u  +  pk)  d  (V) 

V  3 


2  1  3u  V  3v  ^  1  3u  V  3v 

-  ( - +  -  +  — )p  [Peffp  (7 - 7  IF^P  ■"  Pp  '^p]  V 

3  5p  3x  r  3r  ^  'Sp  3x  ^  ' 


2  ,1  3u  V  3v, 

kS3-  _  +_)p  p  V  . 

3  6p  3x  r  3r  ^ 


kS3"  =  kS3  -  (kS3‘)  kp 


=  kS2  +  kS3'‘ 


=  kSl‘  +  kS3‘  . 


-  kS2  +  kS3 


(4-20) 


(4-21) 


(4-22) 

(4-23) 


(4-24) 


Turbulent  Kinetic  Energy  Dissipation  Rate  Equation  Source  Term 


Se  =  —  (Cl  -  C2  p  e)  +  P  e  V*_u 
k 

Su^  and  Spk  are  obtained  by  integrating  the  above  expression  over  the 
control  volume  depicted  in  Figure  4-1. 
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Energy  Equation  Source  Term 
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Integrating  this  expression  over  tne  control  volume  shown  in  figure 
4-1  results  in  Suj  and  Spj. 
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where 


=  volume  of  control  volume  of  previous  time  step 


=  volume  of  control  volume  at  present  time  step. 
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Substituting  the  above  results  into  the  integrated  version  of 


equation  (2-18)  results  in  the  discretized  non-flow  energy  equation: 
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4.2  Special  Procedure  for  Calculating  the  Pressure  Field 


In  order  to  solve  the  u  and  v  momentum  equations  and  the  energy 


equation  the  local  flow  pressure  must  be  specified.  Two  algorithms  exist 


which  lend  themselves  to  calculation  of  P  for  the  present  analysis:  SIMPLE 


[76,78]  and  SIMPLER  [78,87].  Though  the  second  algorithm  has  the  potential 


of  accelerating  the  convergence  rate  of  the  solution  scheme  for  each  time 


step,  it  was  found  arduous  to  use  in  this  analysis.  This  was  due  to  the 


complexity  of  the  source  terms  and  the  time  variations  of  the  boundary 


conditions.  SIMPLE  (Semi  -  Implicit  -  Method  for  Pressure  -  Linked 


Equations)  on  the  other  hand  proved  itself  to  be  "simpler"  to  use  in  this 


analysis  and  was  therefore  incorporated.  For  completeness,  an  outline  of 


the  details  are  given  below  for  the  transformed  flow  field. 
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The  continuity  equation  (equation  4-1)  restated  in  a  less  compact 


form  is: 

no  _ 

Mp  -  Mp  +  (pV)n  Ap  -  (pv)s  Aj  +  {p1J)e  Ag^  -  (pu)y^  Agw  "  0  .  (4-41) 

At  the  outset  of  a  specific  time  iteration,  P,  u  and  v  can  only  be 
guessed.  The  actual  pressure  (P)  at  a  nodal  point  is  equal  to  the  guessed 
pressure  (P*)  plus  the  local  pressure  correction  (P')  and  the  global 
pressure  correction  (P').  The  present  section  deals  with  the  pressure 
correction  due  to  local  velocity  variations  only.  Section  4.3  will  cover 
the  global  pressure  correction.  Therefore,  for  the  present,  we  will 
define  P  as: 


P  =  P*  +  P'  . 


Simi 1 arly. 


u  =  u*  +  u'  , 
U  =  "u*  +  "u' 


and 


V  =  V*  +  V  '  , 


(4-42) 


(4-43) 

(4-44) 


where 


u  =  transformed  velocity  defined  by 
equation  (2-15). 


From  equation  (4-42)  we  can  define  Uw  as  follows: 

ITw  =  Uw*  +  Uw'  , 


and 
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Thus  we  can  express 
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Similarly,  we  can  define; 


Uyy  =  ue*  +  de  (Pp  -  Pe)  , 


vn  = 


+  dn  (Pp  -  Pfj)  , 


Vs  =  V5*  (Ps  "  ^P^  » 


s  ■ 


(4-45) 


(4-46) 


(4-47) 


(4-48) 


where  dw,  de,  dn  and  ds  are  determined  from  the  discretized  momentum 
equations  for  the  control  volume  surrounding  the  appropriate  specified 
flow  field  velocity.  For  example,  the  discretized  equation  for  the  momen¬ 
tum  equation  for  Uyy,  where  Uyy  is  at  the  center  of  the  control  volume  given 
in  figure  4-2,  is  given  by: 
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Vuw  Vg  is  the  volume  of  the  control  volume 

W 

around  Uyy  (see  figure  4-2). 
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Note  that  because  - 


=  guessed  value  of  Uy^  which  is 
labeled  in  figures  4-1  and  4-2. 


and  neither  nor  Up  vary  with  respect  to  a  change  in  the  expression 

(Pw  -  Pp): 

dw  = - =  - 

3(PW  -  Pp)  3(PW  ~  Pp) 

Solving  equation  (4-49)  for  Uyy  and  taking  the  partial  derivative  of 
the  result  with  respect  to  the  expression  (Py  -  Pp)  results  in: 


(4-50) 


aPu  +  -  Sp^ 


Simi 1 ar ly , 


(4-51) 


-  SPu, 


where 


Vug  -  is  the  volume  of  the  control  volume  used  in  the  discretization 

of  the  u-momentum  equation  about  ue  in  an  identical  manner  to  the 
formulation  of  about  uw 

Pug  -  center  nodal  point  of  control  volume  for  ue  which  is  equivalent 
to  the  point  Em  in  figure  4-2. 


For  clarity,  as  shown  in  figure  4-1,  Ue  is  the  velocity  vector 


passing  through  the  east  face  of  the  unstaggered  control  volume  while  Uyy 
passes  through  the  west  face,  Vp  the  north  face  and  V5  the  south  face  of 
this  control  volume.  Identical  formulations  to  the  above  result  in  the 
following  similar  terms  for  Vp  and  V5: 


j'.  V.iVtf 
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Substitution  of  equations  (4-50)  through  (4-53)  into  equations  (4-45) 
through  (4-48)  respectively,  and  substitution  of  the  resulting  four 
expressions  into  equation  (4-41)  yields  in  compact  form  the  following; 
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M  -  subscript  on  Pf/]  which  indicates  capital 

value  of  the  corresponding  little  m  (i,e.,  when 

I  I 

m  =  w.  Pm  =  Py,  etc). 

In  the  overall  solution  scheme,  the  u  and  v  momentum  equations  are 
solved  first  and  all  four  dm's  calculated  for  use  in  equation  (4-54).  Then 
equation  (4-54)  is  solved  and  the  P'  values  obtained  are  used  in  equations 
(4-42),  (4-47),  (4-48)  as  well  as  equation  (4-55)  and  (4-56)  below  to 
obtain  updated  values  of  P,  uw,  we,  Vn  and  vj. 

I  I 

Uw  =  Uw*  +  dw  (Pw  -  pp)  ,  (4-55) 


I  I 

Ug  =  Ue*  +  de  (Pp  -  Pe)  .  (4-56) 

4.3  Perturbation  of  Continuity  Equation 

Inside  internal  combustion  engines,  pressure  variations  are  produced 
by  both  fluid  motion  and  by  the  expansion  and  compression  of  the  air  due 
to  piston  motion.  The  formulation  given  in  section  4.2  describes  the  local 
pressure  correction  due  to  spatial  variations  in  the  fluid  velocities. 

The  procedure  used  to  determine  the  global  variation  of  pressure  due  to 
piston  motion  involves  the  pertubation  of  the  continuity  equation  and  is 
the  subject  of  this  section. 

With  the  assumption  of  a  perfect  gas,  density  is  a  function  of  T  and  P. 

P  =  P  (P.T)  . 

Earlier  researchers  perturbated  the  expression  for  density  with  respect  to 
both  P  and  T  and  obtained  global  correction  terms  for  the  temperature  and 
pressure  using  both  the  energy  and  continuity  equations  [41,67].  These 
formulations  required  a  global  conservation  of  both  energy  and  mass  for 


convergence.  However,  it  has  been  found  that  the  energy  equation,  using 
the  discretization  for  TSl  given  in  equation  (4-32)  will  correct  itself 
globally  for  changes  in  temperature  due  to  piston  motion,  and  that  it  is 
sufficient  to  only  require  a  global  conservation  of  mass  [68].  Thus  the 
expression: 

p  =  p*  +  (i£)  7-  (4-57) 

is  used  in  our  present  analysis  where 

p*  =  guessed  value  of  density  obtained  from  use  of  equation  (2-10), 


P'  =  global  pressure  correction  defined  below  in  equation  (4-58), 

( — )t  =  -  for  a  perfect  gas 

Rgas  T 

and 

n 

Rqas  =  I  fni  Ri  (see  equation  2-10). 

i  =  l 

Note  that: 


8p  V 

=  (P*  +  (— )t  f^’)  — 
3P  at 

=  Mp*  +  0p*  P' 


Substitution  of  this  result  into  equation  (4-10)  yields: 

—  0  .fi  .0 

Mp*  +  0p*  P'  -  Mp  +  mp  -  ms  +  me  -  mw  =  0  . 


Summing  the  above  equation  over  the  entire  flowfield  (F)  and  solving  for 
the  global  pressure  correction  yield: 


This  correction  is  incorporated  into  the  solution  scheme  in  two  places 


It  is  inserted  into  the  source  term  of  equation  (4-54)  so  that  SUp* 


modified  to 

0  *  _ 

SiJp'  =  Mp  -  Mp  -  6p*  P'  -  (pv)n  An  +  (pv)s  Ag 
Also  P"'  is  used  to  correct  P  as  follows: 


(4-59) 

(pu)e  Aew  +  (pu)w  Aew  • 


p  =  p*  +  P' 


(4-60) 


It  should  be  noted  that  the  global  correction  is  made  first  in  the 
iteration  scheme.  Then  equation  (4-54)  is  solved  with  the  equation  (4-59) 
modification  of  Sup,  to  obtain  the  local  pressure  correction  values  (P'). 
Then  equation  (4-42)  is  used  to  adjust  the  values  of  the  flow  field 
pressure  locally. 

4.4  Solution  Procedure 


An  overview  of  the  logic  of  the  code  used  in  the  solution  of  this 
problem  is  given  in  figure  4-4.  The  general  solution  technique  used  to 
solve  the  discretized  form  of  the  governing  equations  for  each  variable  is 
an  iterative  line-by-line  method  in  which  an  initial  guess  of  the  flow 
field  values  is  made  and  then  is  improved  upon  from  one  line  to  another 
[78,88].  The  equations  for  the  points  along  a  particular  line  are 
solved  by  temporarily  assuming  that  the  most  recently  calculated  values  of 
the  terms  on  the  neighboring  lines  are  known,  reducing  the  equation  for 
each  point  along  the  line  (either  N-S  and  E-W)  to  only  three  unknowns. 


Thus  the  equations  along  each  line  take  the  form  of  a  tri-diagonal  matrix 
and  are  solved  using  a  simple  Tri-Diagonal  Matrix  Algorithm 

The  program  which  directs  the  order  of  the  solution  of  the  problem  is 
called  CONTROL.  Upon  running  the  program,  CONTROL  calls  upon  a  series  of 
subroutines  to  set  up  the  initial  geometry  and  specify  initial  property 
values.  Control  then  moves  the  solution  into  a  large  time  loop  which  is 
designed  to  take  the  engine  through  the  desired  number  of  complete  revolu¬ 
tions.  The  program  proceeds  through  this  time  loop  a  time  step  at  a 
time  until  it  reaches  the  final  time  step  (tpinal)-  Within  a  given  time 
step  the  solution  scheme  is  directed  into  a  property  iteration  loop  to 
simultaneously  determine  all  the  properties  at  that  time  using  the  method 
outlined  in  the  previous  paragraph.  The  program  will  not  leave  this  pro¬ 
perty  iteration  loop  until  the  maximum  value  of  the  source  residuals  for 
each  of  the  u,  v,  P',  and  T  equations  is  less  than  .005.  The  values  of 
these  residuals  are  found  by  first  computing  the  difference  between  the 
left  and  right  hand  sides  of  equation  (4-2)  for  all  the  nodal  points  of  a 
given  property  in  the  flow  field.  Then  these  differences  are  summed  over 
the  entire  flow  field  and  divided  by  quantities  appropriate  to  each  parameter 
to  obtain  the  overall  source  residual  for  that  primitive  variable.  Once 
the  largest  of  these  residuals  goes  below  .005  and  last  three  iteration 
values  for  P'  are  all  less  than  .01,  the  program  is  returned  to  the  time 
loop  where  various  old  time  values  are  stored  and  new  geometric  values 
determined  in  preparation  for  the  next  time  step.  In  the  event  the 
program  exceeds  a  fixed  limit  of  iterations  it  is  automatically  advanced 
to  the  next  time  step,  even  though  these  two  convergence  criteria  have  not 
been  met  yet.  At  certain  prespecified  time  steps  along  the  way,  CONTROL 
will  direct  the  storage  of  data  into  files  for  future  analysis  and  plotting. 


Chapter  5 
Valve  Model 

The  problems  encountered  in  the  modeling  of  valve  motion  in  the 
flow  region  of  the  numerical  solution  scheme  were  numerous.  Several 
solution  attempts  were  made  before  the  final  model  contained  herein  was 
arrived  at. 

The  primary  problem  encountered,  that  of  having  a  moving  boundary 
pass  through  solution  nodal  points  which  were  also  moving  due  to  the 
transformation  of  the  governing  equations,  was  solved  by  dividing  the 
flow  field  solution  scheme  into  3  regions  as  mentioned  in  Chapter  2. 

As  can  be  seen  in  Figure  6-1  in  both  region  1  and  region  3,  the 
absolute  distance  of  the  nodal  points  in  the  east-west  direction  from 
their  nearest  north-south  boundaries  do  not  change.  In  region  1  they 
are  slaved  to  the  top  of  the  cylinder  while  in  region  3  they  are  fixed 
to  the  piston  even  though  it  is  moving.  In  region  2  the  geometric  posi¬ 
tions  of  the  nodal  points  in  the  east-west  direction  relative  to  the  top 
of  the  cylinder  are  cha  ging  in  accordance  with  the  transformation  given 
in  section  2.2. 

The  interest  in  developing  a  good  model  for  the  heat  transfer  in 
the  valve  as  well  as  through  the  cylinder  and  piston  linings  resulted  in 
the  development  of  a  method  to  move  a  block  of  nodal  points  representing 
the  valve  back  and  forth  through  the  solution  scheme.  However,  con¬ 
vergence  problems  lead  to  the  partial  adoption  of  the  use  of  an  infini¬ 
tely  thin  valve  reported  by  researchers  in  [69]  instead.  It  was  later 
realized  that  extreme  underrelaxation  procedures  may  have  made  the  use 
of  a  "thick"  valve  in  the  flow  field  feasible.  However,  the  former 


method  was  not  reverted  to,  because,  in  the  interim  analysis,  a  method 
was  discovered  of  giving  the  valve  reasonable  thickness  for  the  purpose 
of  solving  the  conjugate  heat  transfer  problem  while  allowing  it  to  be 
infinitely  thin  for  the  solution  of  the  other  flow  field  parameters. 

The  approach  utilized  herein  therefore  views  the  valve  as  infinitely 
thin  in  relation  to  the  parameters  u,  v,  k,  e,  f,  mf  and  P.  The  valve 
is  viewed  as  having  thickness  from  the  view  point  of  T,  Also,  fluid 
entrainment  due  to  valve  motion  is  considered  in  the  computation  of  the 
unsteady  terms  of  the  flow  field  finite  difference  equations  for  all 
parameters.  Details  of  this  model  are  contained  in  subsequent  subsec¬ 
tions  . 

5.1  Solution  Regions 

Figure  5-1  shows  the  division  of  the  solution  regime  into  four 
solution  regions,  three  flow  field  regions  and  one  non-flow  region. 

Only  in  region  2  is  it  necessary  to  use  equations  in  a  transformed 
state.  In  regions  1,  3,  and  4  the  equations  are  solved  in  an  untrans¬ 
formed  state,  though  special  attention  must  be  given  to  the  area  of 
region  4  adjacent  to  regions  2  and  3  (explained  later). 

In  the  actual  solution  of  the  entire  flow  field,  the  multiple 
region  concept  was  employed  in  such  a  way  as  to  arrive  at  a  simultaneous 
solution  of  the  flow  field  variables  at  all  the  nodal  points  at  one  time 
without  having  to  iterate  back  and  forth  between  the  flow  field  regions. 
However,  iteration  back  and  forth  between  the  non-flow  field  region 
(region  4)  and  the  flow  field  regions  (regions  1-3)  was  required.  This 
simultaneous  flow  field  solution  was  obtained  by  using  the  transformed 
form  of  the  equations  (equation  (2-17))  throughout  all  the  flow  field 
regions  and  by  noting  that  for  regions  1  and  3  the  values  for  the 
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expression  x*Up  used  in  the  computation  of  u  (equation  2-15)  and  in  the 
determination  of  the  transformed  source  term  for  the  energy  equation  Sj 

(see  table  2-1)  are  as  follows: 

region  1:  x*up  =  0;  (5-1) 

region  3:  x*Up  =  Up.  (5.2) 

Use  of  the  transformed  equation  throughout  also  necessitated  the 
use  of  the  absolute  distances  in  the  computation  of  the  areas  (Ap^  and 
and  the  volumes  of  the  control  volumes  as  well  as  all  the  other 
X  coordinates  distances  used  in  the  discretized  terms  of  the  governing 
equations.  In  the  solution  computer  code  the  5p  •  6x  terms  contained  in 
the  governing  equations  were  computed  as  follows  ( 6x  is  any  distance  in 
the  X  coordinated  system,  see  figures  4-1,  2  and  3): 

Regions  1,  2  and  3 

6p  6x  =  6^  ,  (5-3) 


where  for 

Region  1: 

62  =  6x, 

Region  2: 

62  =  6p  4x , 
and 

Region  3; 


(5-3a) 

(5-3b) 


«z  =  (5-3c) 

At  the  boundary  areas  of  the  nodal  points  control  volumes 


corresponding  to  the  interfaces  between  regions  1  and  2  and  regions  2 
and  3  special  attention  had  to  be  given  in  the  application  of  the  trans¬ 
formation  of  the  u-momentum  equations  to  obtain  an  accurate  solution. 


Tnis  was  due  to  the  fact  that  in  the  solution  of  the  finite  difference 
equations  for  the  u-momentum  equation  at  these  interface  nodal  points, 
the  control  volume  boundaries  at  the  interface  can  be  moved  by  our 
transformation  from  one  region  to  another.  See  Appendix  D  for  an 
explanation  of  the  special  procedure  used  in  this  case.  No  special 
problems  arose  in  the  computation  of  the  other  primitive  variables. 

As  mentioned  in  section  2.3  the  solution  of  the  Energy  Equation  in 
region  4  is  carried  out  without  the  use  of  transformed  equations. 

Because  the  portion  of  region  4  adjacent  to  regions  2  and  3  is  changing 
size  and  the  relative  positions  of  the  nodal  points  in  this  portion  of 
region  4  must  be  repositioned.  This  is  necessary  in  order  to  make 
possible  the  matching  of  the  boundary  conditions  at  the  interface  of 
regions  2  and  3  with  region  4.  To  accomplish  this,  a  very  different 
procedure  had  to  be  employed.  See  Appendix  C  for  the  specifics  of  this 
procedure. 

5 . 2  Valve  Inlet  and  Valve  Boundary  Conditions 

To  account  for  the  presence  of  a  functioning  valve,  boundary  con¬ 
ditions  must  be  specified  for  not  only  both  sides  of  the  valve  itself, 
but  also  for  the  orifice  opening  created  when  the  valve  is  open.  The 
thermodynamic  behavior  of  the  gases  due  to  heat  loss  and  heat  absorption 
to  and  from  the  cylinder  and  piston  coupled  with  pressure  and  velocity 
changes  caused  by  piston  motion  slightly  complicate  the  boundary  con¬ 
ditions  at  the  orifice.  However,  with  the  exception  of  the  energy 
equation,  boundary  conditions  at  the  inner  and  outer  surfaces  of  the 
valve  are  easily  specified. 

Orifice  Conditions 

The  nature  of  the  boundary  conditions  at  the  valve  orifice  vary 
depending  on  whether  the  flow  of  gas  is  into  or  out  of  the  engine.  It 
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was  discovered  that  because  of  the  effects  of  heat  transfer  between  the 
gas  and  the  walls,  the  effects  of  pressure  changes  caused  by  combustion 
and  piston  motion,  along  with  other  factors,  the  direction  of  gas  flow 
did  not  always  correspond  with  the  direction  of  piston  motion. 
Obviously,  the  driving  function  is  the  difference  between  the  ambient 
air  pressure  and  the  value  of  the  pressure  at  valve  orifice. 

^atm  ^  Porifice*  orifice  boundary  conditions  are  specified 
as  follows  [69j: 


and 


u  = 


k  = 


2  Cd 


Patm  ■  P 


(.03  u)2  . 

.0075  ryg]yg 


mf  =  0, 

f  =  0, 


T  =  I’atm  (298.  K) , 

where 

Co  =  .7. 

Patm  ^  Porifice  orifice  boundary  conditions  instead  are: 


(5-4) 

(5-5) 

(5-6) 

(5-7) 

(5-8) 

(5-9) 


u  =  - 


2  Cp 

P 


Patm  ■  P 


(5-10) 


3k 

3Z 


0, 


r* 


(5-11) 


v 
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l.'l  V*-' 


.4 


3e 

—  =  0  , 

3z 

(5-12) 

o 

II 

(5-13) 

32 

3f 

—  =  0 

(5-14) 

32 

3T 

-  =  0  . 

32 

(5-15) 

Valve  Surface  Boundary  Conditions 

When  the  valve  is  closed  it  is  acting  in  the  same  geometric  capa¬ 
city  as  the  cylinder  head,  and  with  the  exception  of  the  energy  equation 
the  boundary  conditions  are  those  given  in  section  2.5  for  u,  v,  k,  e, 
mf  and  f  for  the  cylinder  head.  See  equations  (2-24,  25,  26,  40,  41, 

42,  43,  68,  69,  70  and  71). 

When  the  valve  is  open,  for  the  purpose  of  fluid  flow,  it  is  con¬ 
sidered  to  be  an  infinitely  thin  surface.  Boundary  conditions  for  v,  k, 
e,  mf  and  f  on  each  side  of  this  surface  are  specified  in  a  manner  iden¬ 
tical  to  that  used  in  section  2.5  for  both  sides  of  the  valve.  On 
the  interior  side  of  the  valve  they  are  specified  precisely  as  for  the 
case  when  the  valve  is  closed  with  the  exception  that  the  distance  x^g 
(see  3yw+g  in  Figure  2-2)  is  now  determined  from  the  new  valve  surface 

position  to  the  interior  adjacent  nodal  point.  See  Figure  5-2. 

+ 

Obviously  the  near  wall  values  (^ig^^,  pg„,  Xw,  Jw^^  (({>))  used  in 
determing  these  correspond  to  the  properties  of  the  gas  at  the  nodal 
points  immediately  adjacent  to  the  interior  side  of  the  valve.  For  the 
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exterior  side  of  the  valve  (top  of  valve),  the  specification  of  the 

boundary  conditions  for  these  five  parameters  is  identical  to  the  boundary 

conditions  given  in  section  2.5  for  the  piston  surface  with  the 

understanding  that  the  distance  Xgg  and  the  near  wall  property  values 
+ 

(4>ge.  Pge»  ^e«  ‘^W0  properly  determined.  In  other  words,  that 

their  specification  be  based  on  the  fluid  properties  of  the  gas  at  the 
outboard  nodal  points  immediately  adjacent  to  the  external  side  of  the 
valve,  and  the  distance  between  these  nodal  points  and  the  external 
valve  surface. 

^hen  the  valve  is  functioning  the  u-momentum  equation  boundary  con¬ 
dition  for  the  valve  is  specified  as  follows; 

^valve  *  ^valve  (t)  »  (5-16) 

where 

^valve  (^)  =  velocity  of  the  valve  at  time  t  as  specified 
by  the  velocity  profile  for  the  valve  given  in 
Table  5-1. 

Though  the  energy  equation  boundary  conditions  for  the  inter  ior  and 
exterior  sides  of  the  valve  are  also  similar  to  those  used  for  the 
cylinder  head  and  piston  respectively,  they  differ  in  the  manner  in 
which  they  are  inserted.  From  the  view  point  of  the  flow  field,  for  all 
governing  equations,  including  the  energy  equation,  the  valve  is  assumed 
to  be  infinitely  thin,  but  is  considered  to  have  thickness  for  the  pur¬ 
pose  of  solving  the  heat  conduction  equations  within  the  valve.  See 
figure  5-2.  The  procedure  for  accomplishing  this  involves  solving  the 
heat  conduction  equations  for  the  valve  for  a  set  of  special  nodal 


points  located  at  an  imaginary  position  outside  of  the  flow  field.  This 
is  done  in  an  iterative  manner  with  the  energy  equation  for  the  flow 
field.  For  convenience,  the  geometric  position  of  the  nodal  points  used 
to  solve  the  conduction  equations  when  the  valve  is  in  the  flow  field 
(ie,,  is  open)  are  the  same  points  used  as  for  when  the  valve  is  closed. 
In  other  words,  whether  the  valve  is  open  or  closed,  its  nodal  points 
are  considered  to  be  an  extension  of  the  cylinder  head  portion  of  region 
4.  The  thickness  of  the  valve  is  assumed  to  be  69  from  the  viewpoint  of 
the  non-flow  energy-equation  nodal  points  used  to  solve  heat  flow  in  the 
valve  (see  figure  5-2b).  In  the  flow  field,  though  the  valve  is  assumed 
infinitely  thin,  it  is  assumed  to  have  different  temperatures  on  each 
side  of  its  surface  and  The  key  to  accomplishing  the 

simultaneuous  solution  of  the  flow  and  non-flow  equations  is  the  cyclic 
juggling  (during  each  iteration)  of  the  orifice  temperature  (Toriface) 
and  the  temperature  one  nodal  point  to  right  of  the  oriface  (Ti=istep) 
with  the  interior  surface  temperature  of  the  valve  (Ty^p^.)  and  the  tem¬ 
perature  of  the  gas  at  the  nodal  point  immediately  adjacent  to  the 
interior  valve  surface  (^gv-jp^)  respectively.  Thus,  when  the  flow  field 
equations  are  being  solved,  Torifice  and  Ti=istep  are  located  as  shown 
in  figure  5-2a.  Later  in  the  iteration,  when  the  non-flow  equations  are 
being  solved,  the  temperature  values  at  the  nodal  points  corresponding 
to  Toriface  oow  Ty^p^  and  that  corresponding  to  Ti=istep  is  now 
Tgy-nf  Obviously,  the  spacings,  and  other  properties  between  and  at 
these  nodal  points  must  also  be  properly  represented  at  the  appropriate 
times  within  each  iteration.  Rather  than  re-specify  these  spacings  and 
properties  in  the  same  cyclic  manner  as  with  the  temperatures  it  was 
found  convenient  to  determine  and  keep  track  of  the  total  heat  flux 


terms  for  flow  of  heat  from  both  sides  of  the  valve  to  the  flow  field 
gas  (q'‘valvev^  ^nd  q"valvee)-  figure  5.2.  These  heat  flux  terms  are 
used  in  both  the  solution  of  the  flow  and  non-flow  equations  and  are 
updated  once  during  each  iteration.  Following  the  development  and 
nomenclature  pattern  outlined  in  section  2.5,  these  total  flux  terms  are 
specified  below  and  serve  as  the  boundary  conditions  for  the  east  and 
west  surface  of  the  valves  for  both  the  flow  and  non-flow  regions  when 
the  valve  is  open. 


^valvCo  ^  11.63: 


»  I 


t  . 


q  valvee  =  —  (^gvint  '  ’’’vint^- 


•"egv 


(5-17) 


^valvea  ^  11.63: 


q  valveg  “ 
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P  +  P 
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Xvalve^  <  11.63: 


g  valve^^  = 
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^valveyy  ^  11.63: 


q  valvey^  = 
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The  boundary  condition  for  the  equation  along  the  side  (north  border) 


of  the  valve  is  as  follows  when  the  valve  is  open: 


(~)r=valve  radius  "  0  . 

or 


(5-21) 


This  corresponds  to  the  assumption  that  the  end  surface  of  the  valve  is 
adiabatic.  This  assumption  made  it  possible  to  match  boundary  conditions 
for  both  the  flow  region  solution  when  the  valve  is  infinitely  thin  and  the 
non-flow  solution  region  -  where  the  valve  has  thickness.  Also,  this 
assumption  alleviated  the  paradox  of  how  there  can  be  measurable  heat  flux 
through  an  area  which  is  both  zero  and  finite  at  the  same  time. 

In  the  event  the  valve  is  closed  and  thus  is  seated  into  the 
cylinder  head  the  adiabatic  assumption  is  dropped  and  the  boundary  con¬ 
dition  for  the  side  of  the  valve  becomes  that  given  by  equation  (2-65). 

The  valve  exterior  surface  boundary  condition  in  this  case  is  identical 
to  that  given  by  equation  (5-19)  with  the  value  of  used  in  these 

expressions  set  equal  to  atmospheric  temperature. 

5.3  Effects  of  Fluid  Entrainment  Due  to  Valve  Motion  on  Time  Deoendent 


Properties 

In  the  formulation  of  a  finite  difference  scheme  to  accommodate 
valve  motion  through  a  set  of  fixed  nodal  points,  it  was  noted  that, 
extreme  care  was  necessary  in  the  calculating  of  the  values  of  the  prim 
itive  variables  on  either  side  of  the  valve  in  the  current  time  step 
using  the  geometrically  corresponding  values  of  these  properties  at  the 
previous  time  step.  As  the  infinitely  thin  valve  moves  in  one  direc¬ 
tion,  the  values  of  v,  P,  T,  k,  e,  mf,  and  f  along  with  density  located 
immediately  ahead  of  the  valve  "pass  through"  the  valve  and  are  used  as 
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the  old  time  values  in  the  calculation  of  these  variables  on  the  other 


f] 


r* 


side  of  the  valve  in  the  next  time  step.  To  alleviate  this  very 
unrealistic  situation,  the  assumption  was  made  that  the  fluid  property 
values  in  the  boundary  layers  on  both  sides  of  the  valve  are  entrained 
between  time  steps  and  are  available  for  use  in  the  next  time  step  at 
the  next  geometric  nodal  point  in  the  direction  of  valve  motion.  To 
make  this  realistic  the  valve  motion  and  the  nodal  point  spacing  in  the 
vicinity  of  the  valve  had  to  be  closely  coordinated.  To  simplify  com¬ 
putation,  the  valve  was  only  actually  relocated  every  other  time  step  as 
shown  in  Table  6-1,  but  is  considered  to  have  a  velocity  roughly 
corresponding  to  its  rate  of  motion  over  each  two  time  steps.  For  the 
time  step  in  which  the  valve  is  actual ly  moved  geometrically,  the  mathe¬ 
matical  representation  of  the  specification  of  the  old  time  values  for 
use  in  the  current  time  step  are  as  follow: 

For  an  opening  valve  along  both  sides  of  the  valve  surface  - 

4>old  (i.j.h)  =  4i(i-l,  j,  n-1)  ,  (5-22) 

where 

i  =  indicy  used  to  locate  the  position  of  the  nodal  point 
in  the  z  direction, 

i-1  =  indicy  of  next  nodal  point  in  the  solution  scheme 
immediately  to  the  west  of  a  nodal  point  at  i, 

j  =  indicy  used  to  locate  the  position  of  the  nodal  point 
in  the  r  direction, 
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n  =  time  indicy,  indicating  the  current  time  step, 
n-1  =  previous  time  step, 
and 

())  -  understood  here  to  represent  either,  v,  P,  T,  k,  e,  mf 
and  p. 

Similarly,  for  a  closing  valve  the  old  time  values  along  both  sides 
of  the  valve  surface  are  given  by: 

♦old  ( i  =  ♦(’  +  !»  '^-1)  • 


Chapter  6 
FUEL  NOZZLE  MODEl 

In  hypergolic  combustion,  fuel  is  brought  into  the  engine  in  gaseous 
form  under  extremely  high  pressures.  To  date,  in  research  conducted  with 
this  method  of  combustion,  the  same  type  of  orifice  openings  have  been  used 
to  inject  the  fuel  as  are  common  in  diesel  engine  construction  [1,18].  See 
figure  6-1.  At  the  exit  of  the  orifice  in  our  case  the  flow  is  assumed  to 
be  chocked  flow.  This  is  due  to  the  large  pressure  drop  between  just 
inside  the  fuel  injector  and  the  combustion  chamber.  Thus  the  fuel  is 
entering  at  a  Mach  number  equal  to  1.0,  which  presented  numerous  problems 
in  the  modeling  of  the  nozzle.  Rather  than  attempt  to  model  the  actual 
geometric  characteristics  of  the  fuel  injector,  the  object  of  this  analy¬ 
sis  instead  was  to  obtain  a  model  which  gave  a  reasonable  representation 
of  the  fuel  injection  process. 

The  fuel  injector  shown  in  figure  6-l(a)  is  made  up  of  4  orifice 
holes,  each  with  a  diameter  (0)  equal  to  .0105  inches.  The  flow  coming 
out  of  these  orifices  is  such  that  after  a  short  period  of  travel,  the 
streams  from  each  orifice  begin  to  join  and  form  the  shape  of  a  cone  in 
the  manner  shown  in  figure  6-l(b).  To  accommodate  an  axisymmetric  analy¬ 
sis,  it  is  reasonable  to  assume  that  after  a  short  distance  of  travel, 
the  fuel  flow  is  uniform  at  a  given  radius  (i.e.,  invariant  in  0  at  a 
given  radius  where  z,r,0  make  up  the  3-D  coordinates).  Because  of  this, 
it  was  decided  to  model  the  injector  as  an  axisymmetric  convergent  nozzle 
in  a  manner  shown  in  Figure  6-2.  This  was  in  consonance  with  the  recom¬ 
mendations  of  experimental  researchers  in  the  field  [81].  The  dimensions 
of  this  nozzle  are  proportional  to  the  size  of  the  orifices  of  the  injec¬ 
tor  given  in  figure  6-l(b)  in  accordance  with  the  following  relationship: 


(b.)  Fuel  Injector  with  fuel  spray  path 
(ideal  path  for  case  of  no  piston 
or  valve  motion.) 
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Figure  6-1.  Diesel  Fuel  Injector. 


^  /  j" 


(b.)  Nozzle  fuel  injector  with  fuel  spray 
path  (ideal  path  for  case  of  no 
prior  internal  flow  field). 


igure  6-2.  Nozzle  Fuel  Injector  Model. 
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(6-1) 


To  be  compatible  with  experimental  work,  the  stagnation  values  of  the 
temperature  and  pressure  of  the  fuel  entering  the  injector  were  set  at 
2000  psi  and  lOOO'F  respectively.  Also  a  compression  ratio  of  7  to  1  was 
used. 

The  flow  in  the  nozzle  is  assumed  isentropic.  This  assumption  makes 
possible  the  specification  of  the  pressure  at  the  throat  of  the  nozzle 
(see  figure  6-2(a))  [89]. 

Y 

P*  ’  Pi  iTTil  .  (6-2) 


Where 

^Cp 

^  -  lOOO’F. 

v-V 

=  1.0147. 


Also,  for  this  case,  the  throat  boundary  condition  values  for  T,  p,  and  V 
can  be  specified  from  the  following  expressions  [89]: 


and 


T*  =  Ti  [—^1  . 

1 

*  -  r  2  1  T-1 

p*  =  PI  [ - 1 


V*  ARfuel  . 
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Also,  mf,  f,  k  and  e  can  be  specified  as: 


it 

mf 

=  1.0, 

(6-6) 

f* 

=  1.0, 

(6-7) 

k* 

=  0.0, 

(6-8) 

e* 

=  0.0. 

(6-9) 

The  k*  and  e*  values  above  follow  from  the  assumption  of  isentropic 
flow  in  the  nozzle  and  thus  no  turbulence. 

The  viscosity  of  the  fuel  is  estimated  as  [90]: 

ufuel  =  Mfi^^inar  x  10*5  • 

The  problems  created  by  a)  the  geometric  limitations  of  the  nozzle 
(equation  6-1),  b)  the  need  to  specify  all  the  boundary  conditions 
including  the  velocity  at  the  throat  location,  c)  the  need  to  approximate 
the  30*  spray  nozzle,  and  d)  the  Mach  1  velocity  boundary  condition  at 
the  throat  of  the  nozzle,  were  all  eventually  solved  (after  many  unsuc¬ 
cessful  modeling  attempts)  by  the  selection  of  the  model  shown  in  figure 
6-3  subject  to  the  explanation  that  follows.  This  arrangement  enabled  the 
values  of  vnoz  AZpoz  modeled  realistically  without  having  to  be 

majorly  concerned  in  regards  to  grid  spacing  in  the  z  direction  in  the 
vicinity  of  the  nozzle.  The  injector  is  located  far  enough  into  the  flow 
field  so  as  to  be  able  to  deposit  the  bulk  of  the  fuel  in  approximately 
the  same  "location"  as  would  be  accomplished  by  a  spray  nozzle  with  a  30* 
injection  angle  positioned  much  closer  to  the  top  of  the  cylinder.  For 
clarity,  it  should  be  noted  that  after  the  fuel  has  traveled  a  certain 
distance  from  the  injector,  fluid  motion  caused  by  piston  action  becomes 
the  dominant  factor  in  fuel  dispersion  from  that  "location"  on.  Thus,  the 
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fuel  injector  model  concerns  itself  primarily  with  the  preliminary  path 
of  the  fuel  and  in  giving  adequate  specifications  of  all  the  primitive 
variables  at  and  near  the  throat  so  that  they  can  in  turn  be  reasonably 
determined  at  the  deposit  location  at  the  end  of  this  preliminary  path. 
Obviously  a  comparison  of  the  computed  flow  field  with  experimental  data, 
could  be  used  to  "tune"  the  values  of  r^Qz.  z^oz  and  AZnoz» 
model  more  compatible  with  the  actual  physics  of  the  flow  field.  Since 
such  data  does  not  exist,  this  analysis  is  considered  to  be  a  reasonable 
first  cut.  The  problem  of  co-specifying  boundary  conditions  for  v,  along 
with  k,  e,  T,  mf,  f  and  P  al 1  at  and  near  the  throat  of  the  nozzle  was 
overcome  by  varying  the  discretization  of  the  basic  transport  equation  at 
the  control  volume  adjacent  to  the  nozzle  as  explained  in  the  paragraphs 
that  follow.  The  problem  with  a  velocity  of  Mach  1  was  solved  by 
modifying  several  of  the  expressions  to  make  them  physically  compatible 
with  the  unique  flow  conditions  in  the  vicinity  of  the  nozzle.  This  also 
will  be  explained  further  in  later  paragraphs  in  this  chapter. 

In  the  finite  difference  solution  of  the  energy  equation,  infor¬ 
mation  at  nodal  point  (1)  shown  in  figure  6-3  is  not  passed  to  nodal 
point  (2),  because  the  boundary  condition  of  the  control  volume  boundary 
corresponding  to  the  nozzle  throat  (*)  is  known.  (Obviously,  due  to  phy¬ 
sics  of  flow,  we  would  not  want  to  use  point  (1)  in  the  solution  of  (2) 
anyway).  These  two  nodal  points  are  decoupled  by  setting  the  terms  for 
the  total  flux  between  them  (aj^  given  in  equation  4-2e)  equal  to  zero. 

A  correct  property  flux  term  which  reflect  the  throat  (*)  total  flux  value 
is  inserted  into  the  solution  of  the  governing  equation  at  nodal  point 
(2)  through  the  use  of  the  source  terms  Su^  and  Spij,  [76].  Specific 
changes  to  the  discretized  equation  contained  in  4.1.2  for  use  in  solution 


of  temperature  at  nodal  point  (2)  are  summarized  in  Table  6-1. 


The  values  of  Uy^  and  ug  corresponding  to  nodal  point  (2)  in  figure 
6-3  (also  see  figure  6-4)  are  set  equal  to  zero.  This  is  done  for  two 
reasons.  First  of  all,  it  is  necessary  to  set  both  values  equal  to  zero 
in  order  to  ensure  that  the  spray  exiting  the  nozzle  is  actually  perpen¬ 
dicular  to  the  nozzle  assembly  as  desired.  Secondly,  a  flow  pattern  simi¬ 
lar  to  that  which  results  from  the  vena  contracta  effect  [76]  is  assumed 
for  flow  between  the  throat  and  point  (2).  For  such  a  flow,  there  is 
assumed  to  be  no  flow  across  the  outer  border  of  the  "modified"  vena 
contracta  boundary  (see  figure  6-4),  and  a  constant  mass  flow  rate  (m) 
within  its  borders.  Though  this  effect  does  not  actually  dictate  u^  and 
Ug  are  both  zero,  it  dictates  that  the  normal  velocity  on  the  borders  of 
the  vena  contracta  must  be  zero.  Specifying  Ug  and  u^  velocities  for 
nodal  point  (2)  to  be  zero  not  only  approximates  this,  but  also  ensures 

that  the  mass  flow  rate  from  the  throat  (*)  to  point  (2)  within  the 

.n 

modified  vena  contracta  are  equal.  Additionally,  it  ensures  that  m^^^ 

(the  mass  rate  of  flow  through  the  southern  boundary  of  the  control  volume 

for  the  solution  of  V3  (control  volume  shown  in  small  dashed  lines  in 

figure  6-4)  is  the  same  as  the  mass  flow  rate  of  fuel  at  the  throat  (*). 

Note,  that  due  to  equal  spacing,  this  area  (Ac,,,  )  passes  through  point 

v(  3) 

(2).  This  arrangement  enables  os  to  modify  the  discretization  for  the  v- 
momentum  equation  in  the  control  volume  for  v(3)  in  a  manner  identical  to 
the  modification  made  to  the  energy  equation  at  point  (2).  This  arrange¬ 
ment  also  enables  us  to  specify  the  values  of  mf  and  f  at  point  (2), 
namely  - 


and 


(6-11) 


1 


i 


f(2)  =  1.0. 

For  the  case  of  chocked  flow,  which  is  assumed  here,  the  expansion 
process  which  takes  place  as  the  fuel  leaves  the  nozzle  is  not  isentropic. 
Our  k-e  equations  are  not  adequate  for  the  computation  of  k  or  e  in 
the  vicinity  of  nodal  point  (2).  It  was  determined  best  to  assume  the 
flow  to  have  reached  a  reasonable  percentage  of  its  full  turbulence  value 
at  point  (2)  and  to  specify  appropriate  values  for  k  and  e  at  this  point. 
Table  5-1  summarizes  the  specific  changes  made  to  the  discretization  of 
the  equations  used  in  the  solution  of  the  primitive  variables  at  nodal 
point  (2)  and  its  vicinity. 

As  can  be  seen  from  this  table,  the  values  of  Su^  and  Spij,  are 
modified  to  incorporate  these  changes.  The  effects  of  compressibility  are 
greater  in  the  control  volume  adjacent  to  the  nozzle.  Thus,  the  values 
Sup'  and  Spp-  are  modified  as  shown  in  Table  6-1  to  give  more  accurate 
treatment  along  the  northern  area  of  the  control  volume  of  nodal  point  (2) 
in  the  equation  scheme  used  to  compute  the  local  pressure  correction  (P'). 
The  expressions  given  in  Table  6-1  for  P'  follow  from  improving  the 
specification  of  density  at  this  border  area  to: 

★  ' 

Pn  =  Pn  +  Pn  .  (6-12) 

where  p*  is  the  guessed  density  and  p'  is  a  density  correction.  Note  that 

p'  can  be  expressed  as: 

LL 

RT  ’ 
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djacent  to  the  throat  opening  of  the  fuel  injector  shown  in  Figure  6-4 


nts  to  the 
t  [1]  ire  se*. 


CHAPTER  7 


ERROR  ANALYSIS  OF  THE  SOLUTION  SCHEME 


In  the  actual  descriti zation  of  the  complex  transformed  governing 
equations  used  in  this  analysis  there  is  much  chance  for  error.  There  is 
also  the  possibility  that  mistakes  can  be  made  in  the  programming  of  the 
solution  code.  Additionally,  even  if  these  first  two  evolutions  are  carried 
out  flawlessly  there  is  the  question  of  just  how  accurately  the  discretized 
equations  can  be  expected  to  solve  the  problem.  The  present  chapter  deals 
with  these  concerns  by  applying  a  novel  and  very  powerful  error  analysis 
technique  developed  by  Shih®^. 

The  analysis  was  conducted  in  such  a  way  that  the  accuracy  of  the 
governing  equations  could  be  tested  both  one  at  a  time  and  in  combination 
with  any  number  of  the  other  primitive  variable  equations  for  the  non¬ 
combustion  case.  It  was  conducted  for  both  the  untransformed  and  transformeo 
versions  of  the  governing  equations  and  represents  the  first  known  error 
analysis  of  this  type  to  the  power  law  descritization  formulation  of  Patankar 
used  herein.  [76] 

7.1.  Error  Analysis  Technique 

The  discretization  of  the  general  transformed  transport  equation  for  the 
primitive  variable  <p  is  given  by  equation  (4-2)  in  section  4.1.2.  This  repre¬ 
sents  an  algebraic  approximation  of  the  actual  governing  equation  for  ^  given 
in  equation  (2-17)  and  can  be  use  to  give  a  good  approximation  to  the  actual 
solution  of  (2-17).  In  real  life,  there  is  no  analytic  solution  to  equation 
(2-17)  for  any  of  its  41's.  However,  suppose  we  specify  that  for  the  flow 
field  equations  in  general  - 


<t>  =  4>  (x,r,t) 


(7-1) 
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subject  to  Dirichlet  boundary  conditions  for  all  primitive  variables  41. 
Obviously  4  (x,r,t)  can  not  satisfy  equation  (4-2)  or  equation  (2-17). 
But  it  can  be  made  to  solve  the  modified  version  of  equation  (4-2)  given 
below; 


no  n 

'"P,  *  %  -  Sp»> 


n  00 

I  +  Mp^  4p^  +  C.Ranal^ 

m 


(7-2) 


where 


and 


f' 

1  when  error  analysis  utilized 
in  computer  code  subject  to 
Dirichlet  boundary  conditions. 

I 

0  For  solution  of  real  problem 
subject  to  boundary  conditions 
specified  in  Chapter  2.5. 


^anal,  =  defined  by  the  below  development 


First  we  define  V4  as  follows: 


_  1 


1 


6p  3t 


9 „ 


=  4  —  (p  6p«)  +  Lpu4  '  6^-^  ]  +717  LrpV4  -  r  r,. 


Substitution  of  this  expression  into  equation  (2-17)  reduced  the 
transformed  version  of  the  general  transport  equation  into  the  following 
compact  expression; 

V4  =  (7-3) 

If  we  boldly  substitute  the  analytic  expression  for  4,  namely 


4  =  4(x,r,t)  into  equation  (7-3),  we  will  obviously  not  obtain  a  solution. 


because  (7-3)  has  no  analytic  solution.  Proceeding  anyway  we  find  upon 
substitution  that 

V(^(x,r,t)  =  (x,r,t)  ^  f^anaU  ,  (7-4) 

where 

Ranal^  =  7$(x,r.t)  -  S<^(x,r,t)  • 


Thus,  i))(x,r,t)  is  the  analytic  solution  of  new  equation  (7-4). 
Discretization  of  equation  (7-4)  using  precisely  the  saTie  technique  as  was 
employed  in  section  4.1.2  results  in  equation  (7-2)  with  C=l.  If  we  spe¬ 
cify  Dirichlet  boundary  conditions  we  can  use  equation  (7-2)  to  numeri¬ 
cally  solve  equation  (7-4)  and  can  then  compare  the  results  of  the 
computation  with  the  "known"  (analytic)  solution  at  each  of  the  nodal 
points  in  the  solution  scheme  See  figure  7-1  for  a  specification  of 
these  boundary  conditions . 

7.2  Formulation  of  Error  Analysis  Test 

After,  much  trial  and  error  and  significant  analytic  analysis,  it  was 
elected  to  make  the  following  specif ications  for  for  the  non-combusting 
case: 

u  s  tx^r  ,  (7-5) 

V  =  -2/3  txr2  ,  (7-6) 

k  =  C6  x*r  (x  +  r)t  ,  (7-7) 

e  i  cg^  x^  r^t  C7  ,  (7-8) 

and 

T  =  t  x2  r2  .  (7-9) 

where 


and 


C5  =  arbitrary  constant  (value  of  .75  used) 


C7  =  arbitrary  constant  (value  of  .8888  used). 


I  ‘  'i 

'  V  ■  ■  *1 


/  / 


The  selection  of  p  =  1  was  also  made  corresponding  to  incompressible 

K 

flow.  The  value  for  R  used  was  that  of  the  gas  constant  for  air.  It 
should  be  noted  that  the  selection  of  equation  (7-5)  in  conjunction  with 
equation  (7-6)  was  necessary  in  order  to  satisfy  the  unsteady  transformed 
continuity  equation  (equation  2-16).  It  was  found  that  the  above  for¬ 
mulation  coupled  with  the  nature  of  the  SIMPLE  algorithm  resulted  in  a 

aP 

decoupling  of  the  primitive  variable  P  from  the  solution  scheme.  The" 

dX 

term  in  the  u-momentum  equation  for  example  was  retained  as  a  negative 
expression  in  Su^j  and  as  a  positive  expression  in  Ranaly  3rid  w^s  thus 
cancelled  out. 

The  substitution  of  equations  (7-5)  through  (7-9)  into  equation  (7-3) 
resulted  in  the  expressions  for  Ranal^’  ^^anal^,  Ranalk’  l^anale’ 

Ranalj  summarized  in  Table  7-1. 

For  comparison  with  Ranalj  various  4>'s  the  expression  Rpum^  'vas  derived 

from  equation  (4-2)  with  the  following  result: 

1  n  0  n 

Rnum  =  —  [ap^  +  Mp^  -  Sp^) 


-  (I  +  Mp^  <^PJ]  • 


(7-10) 


Results  of  the  error  test  conducted  on  the  solution  code  are  discussed 
in  Chapter  8. 
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CHAPTER  8 


DISCUSSION  OF  RESULTS 

In  its  development,  the  solution  code  was  configured  to  handle 
five  specific  engine  cases:  (1)  the  open  orifice  case,  (2)  the  con- 
I  tinuously  closed  valve  case,  (3)  the  open  fixed  valve  case,  (4)  the 

fully  operating  valve  case,  and  (5)  the  fully  operating  valve  with  com- 
bustioh  case.  The  present  chapter  will  report  on  the  results  obtained 
,  from  the  code  for  cases  (3)  and  (5)  only  along  with  results  obtained 

I 

from  the  application  of  the  error  analysis  formulation  of  Chapter  7  to 
the  code. 

8 . 1  Error  Analysis  Results 

Several  series  of  error  test  runs  were  conducted  during  the  develop 
ment  of  the  solution  code  as  a  means  of  checking  for  mistakes  both  in 
the  discretization  of  the  governing  equations  and  in  the  actual 
programming  of  the  code.  The  first  set  of  tests  was  made  on  the  code 
representation  of  the  untransformed  equations  for  u,  v,  k,  e  and  T.  In 
these  tests  the  spacing  was  uniform  for  both  the  z  and  r  coordinate 
systems.  Afterwards,  the  code  representation  of  the  transformed 
equations  for  these  same  five  primitive  variables  was  tested  with  non- 
uniform  spacing  in  the  z  direction  but  with  uniform  spacing  in  the  r 
direction.  The  error  tests  reported  on  in  this  section  are  for  a  series 
of  runs  made  much  later,  in  the  weeks  just  prior  to  the  final  data  runs 
for  engine  cases  3  and  5.  The  rationale  for  conducting  this  final  round 
of  tests  was  to  ensure  that  no  bugs  had  been  inadvertantly  inserted 
into  the  program  in  the  interim. 

The  geometric  values  used  in  the  formulation  of  the  boundary  con¬ 
ditions  for  these  runs  (see  figure  7-1)  are  as  follows: 


'"smal  1  “  1-0144  m, 
■"large  =  1.06424  m. 


distl  =  .04  m 
and 

dist2  =  .0033333  m. 

Values  used  in  the  solutions  of  equations  (B-1)  and  (B-2)  for  Zp 
and  Up  respectively  are  as  follows; 

t  =  .01  m, 

L  =  .06  m, 

C  =  1.0  m  (Clearance) 

and 

RPM  =  1900. 

The  values  of  6p  (obtained  from  Zp  -  see  equation  (2-12))  and  up  define 
the  transformation  [see  equations  (2-13)  through  (2-15)]. 

The  time  domain  boundary  conditions  corresponding  to  the  previous 
time  step  were  specified  for  (|i(x,r,t  -  At)  using  a  time  step  of 

1 

4t  =  Wm  • 

The  value  of  time  used  in  the  analysis  was 


3000 

RPM 


1.57895  sec  . 


As  mentioned  in  Chapter  7  the  error  analysis  scheme  was  programmed 
in  such  a  way  as  to  facilitate  the  testing  of  the  code  representation  of 
one  primitive  variable  at  a  time,  as  well  as  any  combination  of  primitive 
variables  at  a  time.  Obviously,  those  variables  not  being  tested  had  to 
be  specified  analytically  due  to  the  fact  that  the  governing  equations 
are  coupled. 
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In  Appendix  E,  results  of  seven  of  these  runs  are  listed.  Tables 
E-1  through  E-5  represent  individual  variable  runs  for  the  primitive 
variables  u,  v,  k,  e  and  T  respectively.  Table  E-6  represents  a  com¬ 
bined  run  for  all  five  of  these  variables.  For  the  runs  depicted  in 
Tables  E-1  through  E-6  the  r  spacing  is  uniform  while  the  z  spacing  is 
non-uniform.  The  last  run,  given  in  Table  E-7,  was  for  the  combination 
of  primitive  variables  u,  v  and  k.  During  this  run,  the  spacing  in  both 
the  r  and  z  directions  is  non-uniform.  As  can  be  seen  from  data  con¬ 
tained  in  these  seven  tables  the  comparison  between  the  numerically  com¬ 
puted  values  of  the  various  (ji's  with  the  actual  (analytic)  values  is 
excellent  for  all  five  primitive  variables  both  when  they  are  computed 
individually  and  in  combination. 

8.2  Case  3  Results 

A  computation  was  conducted  for  a  dense  mesh  of  23  x  29  for  the  up 
stroke  of  a  two  stroke  engine  with  its  valve  in  a  fixed  position  (open 
orifice).  The  purpose  of  the  computation  was  to  compare  the  ability  of 
the  code  to  compute  the  u-component  of  velocity  with  the  code  of  another 

C  Q 

researcher,  Ramos  whose  valve  model  is  similar  to  the  model  used 
herein.  Ramos'  work  for  the  case  of  a  2-stroke  engine  gave  fair 
agreement  with  recent  experimental  data.  The  objective  of  the  run  was 
the  comparison  of  the  u-velocity  profile  along  a  specific  z  position 
(z  =  .003350  m)  between  the  valve  and  the  top  of  the  cylinder. 

The  engine  along  with  the  valve  model  was  configured  in  a  manner 
identical  to  that  used  in  reference  [69].  The  mesh  size  and  nodal  point 
spacing  of  the  two  computations  were  similar  (22  x  30  mesh  used  in 
[69]).  In  Ramos'  calculation,  the  engine  was  run  through  several 
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complete  cycles  until  periodic  solutions  of  the  flow  field  variables 
were  obtained.  Then  using  the  final  values  of  the  previous  cycle  as 
initial  conditions,  he  made  his  data  computations.  His  work  utilized  a 
UNIVAC  1100/40  and  a  CDC  6600,  and  required  35  and  10  minutes  of  CPU 
time  respectively.  The  present  computations  were  made  on  a  VAX  11780 
and  required  approximately  two  days  of  on  line  time  to  complete  the 
computations  for  a  single  stroke  for  the  23  x  29  mesh  size.  Because, 
obtaining  warm  up  conditions  for  this  calculation  was  not  possible,  it 
was  decided  best  to  compare  only  the  up  stroke  calculations  for  the  u- 
velocity  at  the  z  location  approximately  halfway  between  the  fixed  valve 
and  the  top  of  the  cylinder.  The  comparison  was  made  at  a  crankshaft 
angle  of  300°.  Part  of  the  rationale  for  this  choice  of  comparison  was 
the  expectation  that  the  valve  would  dominate  the  flow  in  this  region 
and  the  lack  of  start-up  conditions  would  have  little  effect  on  the  com¬ 
putation  of  the  u-velocity  at  this  axial  location. 

The  following  are  characteristics  of  the  engine  configuration  used 
in  this  run: 

i  =  .0381  m  (semi-stroke  length), 

L  =  .2032  m  (connecting  rod  length), 
valve  radius  =  .00785  m, 
cylinder  radius  =  .03835  m, 

C  =  .0127  m  (clearance). 


and 


Zy  =  .00793  m  (valve  position). 


Additional  values  included: 

a  =  300°  (crank  shaft  angle). 
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The  run  was  started  at  0  =  180*  with  the  values  u,  v,  k  and  c  set 
equal  to  zero.  The  air  temperature  was  specified  at  300  K  throughout 
and  the  initial  pressure  set  at  atmospheric  pressure.  The  wall  tem¬ 
perature  boundary  conditions  were  set  at  350  K  for  the  cylinder  head, 
cylinder  side  and  piston  surface. 

There  are  only  several  real  differences  in  the  solution  formulation 
of  reference  [69]  and  the  present  analysis.  Ramos'  code,  based  on  the 
Teach-T  code  [78],  used  a  hybrid  central/upwind  spacial  differencing 
scheme  in  computing  the  total  flux  terms  while  as  explained  in  Chapter 
4,  the  present  code  uses  a  power  law  spacial  differencing  scheme  [76]. 

The  present  work  specifies  the  turbulent  kinetic  energy  boundary 
conditions  in  a  manner  identical  to  Teach-T  while  Ramos  assumed  the 
flux  of  zero  for  this  primitive  variable  between  the  wall  and  the  near 
wall  nodal  points.  Also,  the  present  run  used  a  fixed  value  for  the 
specific  heat  throughout.  In  Appendix  E,  figure  E-1  gives  a  vector  plot 
of  the  velocities  computed  in  the  case  3  run.  A  recirculation  zone  is 
visible  behind  the  valve.  Figure  E-2  shows  a  comparison  of  the  u-velocity 
profiles  for  the  z  =  .003350  positions  of  the  two  runs.  For  completeness, 
profile  comparisons  are  given  near  the  z  =  .017037  and  at  the  z  =  .026695 
positions  also.  As  can  be  seen,  there  is  excellent  comparison  for  the 
two  profiles  at  the  z  =  .003350  position.  The  profiles  at  the  other 
positions  do  not  compare  well,  as  expected.  This  is  due  to  the  fact 
that  during  the  intake  stroke,  a  strong  recirculation  zone  is  formed 


below  the  valve.  During  the  up-stroke,  the  effects  of  this  zone  are 
still  present  causing  the  irregular  patterns  shown  for  Ramos'  profiles 
near  the  z  =  .017037  and  at  the  .026695  positions.  Because  the  present 
analysis  assumed  the  velocities  to  be  zero,  throughout  the  flow  field  at 
the  9  =  180°  position  the  effects  of  this  recirculation  are  not  present. 
Thus  the  profiles  at  these  z  locations  are  much  flatter. 

Figure  E-3  shows  the  isotherms  for  this  computation.  As  can  be 
seen  the  temperature  gradient  near  the  walls  are  much  greater. 

8.3  Case  5  Results 

Two  computations  (runs  1  and  2)  were  conducted  for  the  compression 
stroke  and  fuel  ignition  portion  of  the  strokes  to  test  the  nozzle 
model.  One  computation  of  the  first  56°  of  the  exhaust  stroke  (run  3) 
was  conducted  to  test  the  valve  model.  The  computational  mesh  for  Run 
number  1  was  20  x  27  from  the  viewpoint  of  the  energy  equation 
(including  both  the  blow  and  non-flow  regins)  and  17  x  21  from  the  view¬ 
point  of  the  other  flow  field  equations.  The  engine  geometry  used  in 
this  run  is  identical  to  that  utilized  in  case  3  with  the  exception  of 
the  use  of  a  valve  radius  of  .0086  m  and  the  use  of  an  RPM  of  1900. 

Note  that  the  corresponding  time  step  for  this  RPM  was  .0005263  seconds. 
The  engine  was  started  at  a  crank  case  angle  of  180°  corresponding  to 
the  beginning  of  the  compression  stroke.  The  values  of  u,  v,  k  and  e 
were  initialized  at  zero.  The  temperature  and  pressure  specified 
throughout  were  standard  atmospheric  values.  The  temperature  of  the 
valve  was  set  at  400  K  throughout.  Constant  values  of  specific  heat 
were  used  during  the  compression  stroke,  and  then  computed  as  per 
equation  (3-10)  during  combustion.  kSl’  was  set  at  zero  at  the  symmetry 
axis  boundary  nodal  points  (equation  4-18).  Just  as  the  piston  reached 
TDC  (9  =  360°),  the  fuel  injector  was  instantly  inserted  into  the  engine 


and  sprayed  in  fuel  in  a  gaseous  state  (T  =  805  K,  P  =  821,000  N/m  )  for 
5  time  steps.  During  the  compression  stroke  the  mass  fractions  of  the 
working  substance  were  fixed  at  their  atmospheric  values,  and  the 
equations  for  f  and  mf  were  not  solved.  At  and  after  0  =  360°  these 
additional  equations  are  utilized  for  the  power  stroke  calculations. 

Run  2  was  identical  to  run  1  with  several  exceptions.  The  mesh 
size  used  was  17  x  20  and  20  x  26  for  the  flow  field  and  extended  energy 
equation  solution  realms.  The  value  of  kSl'  was  specified  at  the  sym¬ 
metry  axis  nodal  points  as  per  equation  (4-18)  during  compression  and 
set  at  zero  during  combustion.  The  values  of  specific  heats  were  spa- 
cially  computed  during  the  entire  run.  The  temperature  boundary  con¬ 
ditions  used  were  as  specified  in  figure  (2-4).  For  run  1  the  7  nodal 
points  corresponding  to  the  upper  left  hand  corner  of  the  cylinder  con¬ 
figuration  given  in  figure  (2-4)  had  values  of  325°  during  much  of  the 
calculation  instead  of  the  values  shown.  Fuel  was  injected  for  12 
degrees  (two  time  steps).  The  value  of  Prt  for  run  2  was  set  at  1.0 
while  for  run  1  it  was  specified  at  0.9.  Also  the  valve  radius  was 
changed  to  7.5  cm.  For  the  exhaust  stroke  calculations  (run  3),  the 
dependent  variables  were  initialized  identically  to  run  2  and  the  same 
mesh  spacing  and  engine  dimensions  were  used.  Piston  motion  was  begun 
at  0  =  534°  and  continued  for  10  time  steps.  The  valve  was  operated  as 
per  equation  (5-16).  Results  for  these  3  runs  are  given  in  Appendix  E. 
Values  displayed  in  the  Appendix  are  in  SI  units. 

Run  1  results  are  given  in  figures  E-5,6,7,11,12,13  and  14.  Prior 
to  combustion  the  temperature  of  the  flow  field  is  essentially  uniform 
as  shown  in  figure  E-4  and  the  flow  is  in  the  axial  direction  only 
(figure  E-11).  The  duration  of  the  fuel  injection  for  this  run  was  30° 
of  crankshaft  angle  movement.  The  iteration  limit  per  time  step  was 
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specified  as  650  which  was  sufficient  for  all  but  the  time  step 
corresponding  to  6  =  360°.  Figure  E-5  through  E-7  show  in  succession, 
the  growth  and  expansion  of  high  temperature  regions  to  the  right  and 
left  of  the  path  of  the  fuel  spray  for  e  =  360,  372  and  384  respec¬ 
tively.  These  figures  are  given  in  a  modified  x  coordinate  system  in 
order  to  show  the  details  of  the  profiles  in  the  expanding  and 
contracting  portion  of  the  solution  scheme  more  clearly.  The 
corresponding  velocity  vector  plots  for  these  values  of  0  are  given  in 
figure  E-12  through  14.  As  can  be  seen  in  figure  E-12,  there  is  a 
fairly  strong  flow  along  the  side  of  the  cylinder  in  the  direction  of 
the  piston  even  though  the  piston  velocity  is  zero.  Also,  two  recir¬ 
culation  zones  are  created  above  and  below  the  path  of  the  fuel  spray. 

Run  2  results  are  shown  in  figures  E-4,8,9,  15,  16  and  18-25.  For 
this  run,  the  iteration  limit  per  time  step  was  increased  and  the 
source  term  convergence  criteria  was  obtained  for  all  values  of  0. 

Though  the  value  of  P'  did  not  go  below  .01  for  0  =  360°,  it  attained 
acceptably  low  levels.  The  additional  iterations  allowed,  account  for 
the  difference  in  the  isotherm  plot  for  this  run  at  0  =  360°  (figure 
E-8)  and  the  isotherm  plot  for  run  1  at  the  same  crankshaft  angle 
(figure  E-5).  In  figure  E-16  one  can  see  that,  even  though  the  piston 
velocity  is  small,  there  is  a  relatively  large  velocity  component  at  the 
side  of  the  cylinder  wall  in  the  direction  of  the  piston.  This  is  due 
largely  to  the  influence  of  the  fuel  stream  as  well  as  the  fact  that  the 
off  centeredness  of  the  fuel  injector  in  the  flow  fiel  at  0  =  360° 
coupled  with  the  upward  momentum  of  the  flow  due  to  the  piston  motion  of 
the  compression  stroke  caused  the  flow  to  turn  downward  in  the  upper 
corner  of  the  side  wall  at  that  time  (see  figure  E-15).  This  high  velo¬ 
city  value  near  the  wall  accounts  in  part  for  the  high  values  of  the 
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6  =  366°  as  shown  in  figures  E-21  and  E- 23  respectively.  In  figure 
E-20  an  isobar  plot  of  the  specific  heats  for  the  flow  field  at  e  =  366° 
is  given.  As  shown,  the  highest  values  are  located  in  the  vicinity  of 
high  fuel  mass  fractions  (see  figure  E-19)  and  higher  flow  field  tem¬ 
peratures.  The  mass  fractions  of  the  products  at  e  =  366°  shown  in 
figure  E-18  are  high  in  the  vicinity  of  the  recirculation  areas  depicted 
in  figure  E-16  and  conform  also  to  the  high  temperature  regions  shown  in 
figure  E-9.  The  maximum  values  of  both  turbulent  kinetic  energy  and 
turbulent  energy  dissipation  rate  for  6  =  366°  are  located  in  the  vici¬ 
nity  of  the  fuel  injector  as  expected  (see  figures  E-21  and  E-23) .  This 
is  due  primarily  to  the  high  values  specified  at  the  nodal  point  adja¬ 
cent  to  the  injector  orifice  opening  in  the  modeling  of  the  fuel  injec¬ 
tor.  Another  factor  accounting  for  the  high  values  of  turbulent  kinetic 
energy  at  the  side  wall  of  the  cylinder  shown  in  figure  E-21  may  be  the 
inadequacy  of  the  wall  boundary  condition  assumption  made  in  Chapter  2 
in  the  discretization  of  the  k  source  term  (see  equation  2-38).  Plots 
of  pressure  versus  crankshaft  angle  and  the  rate  of  heat  transfer  bet¬ 
ween  the  surface  of  the  piston  and  the  engine  gases  per  crankshaft  angle 
are  given  in  figures  E-24  and  E-25. 

The  results  of  run  3  are  given  in  figures  E-10  and  17.  The 
influence  of  the  valve  on  the  temperature  of  the  flow  field  is  obvious 
in  figure  E-10.  The  velocity  of  the  flow  field  is  essentially  axial 
between  the  piston  and  the  area  below  the  the  valve  (see  figure  E-17) 
and  is  dominated  in  the  upper  portion  of  the  engine  cylinder  by  valve 
motion  and  the  valve  orifice. 

8.4  Concluding  Remarks 

The  model  described  in  Chapters  1  through  6  brings  to  a  focus  on  a 
potentially  significant  engine  configuration  an  assortment  of  recent 


developments  in  the  areas  of  fluid  mechanics,  combustion  and  heat 
transfer.  The  results  outlined  in  sections  8.1  through  8.3  demonstrate 
the  accuracy  and  utility  of  this  model.  It  is  now  a  viable  tool  for  use 
in  the  study  of  low  heat  rejection  engines.  Though  in  its  present  con¬ 
figuration,  it  is  designed  specifically  for  the  study  of  heat  transfer 
in  an  engine  with  a  weak  heat  barrier  in  the  piston  and  using  hypergolic 
combustion,  it  was  developed  in  such  a  way  as  to  facilitate  its  use  in 
other  engine  configurations.  Engines  with  geometries  of  greater 
complexity,  with  different  combustion  conditions,  with  barriers  in  other 
locations  in  the  engine,  etc.  could  be  solved  by  it  without  major  modifi¬ 
cations.  In  the  future  the  code  may  prove  valuable  for  cross  comparison 
purposes  in  the  evaluation  of  other  codes. 

The  code  represents  a  compromise  between  accuracy  and  computational 
efficiency.  In  order  to  obtain  more  precise  wall  heat  transfer  calcula¬ 
tions,  a  low  Reynolds  number  model  with  a  refined  mesh  (such  as  Jones  and 
Launder's  in  [93]  is  needed  for  the  near  wall  calculations.  Real  engines 
are  not  axisymmetric  and  require  three-dimensional  calculations  which 
include  swirl  to  simulate  them.  Though  hypergolic  combustion  shows  pro¬ 
mise  in  future  engine  applications,  most  engines  use  other  combustion 
techniques  which  require  more  sophisticated  combustion  models  to  ade¬ 
quately  simulate.  Also,  a  more  accurate  specification  of  the  specific 
heat,  more  sophisticated  equations  of  state,  higher  order  closure  models 
for  turbulence,  etc,  could  be  incorporated  -  all  increasing  the  required 
computational  time. 

The  code  assumes  the  adequacy  of  the  k-e  model  in  the  rapid 
compression  and  expansion  situations  occurring  in  IC  engines  which  is  not 
the  case.  Turbulence  in  the  engine  flow  field  is  not  isotropic  as  the  k-c 
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model  assumes.  Reference  [94]  points  to  the  need  for  modifying  the  k-c 
model  under  these  circumstances  and  postulates  such  a  modification.  As 
conclusive  improve-ments  are  made  in  the  k-e  model  for  engine  applications 
the  code  presented  herein  should  be  updated.  Also,  as  alluded  to  in  sec¬ 
tion  8.3,  further  study  of  the  boundary  conditions  for  k  at  the  near  wall 
nodal  point  is  in  order. 

In  its  present  state,  the  code  should  hopefully  serve  as  a  baseline 
for  the  analysis  of  the  engine  processes  involved  in  internal  regenera¬ 
tion  and  be  of  aid  in  the  development  of  the  heat  barrier  piston  engine. 

It  may  also  be  of  assistance  in  the  evaluation  and  design  of  engines 
using  hypergolic  combustion. 

Follow  on  research  should  be  conducted  to  make  the  model 
usable  in  the  study  of  stronger  heat  barriers.  In  such  conditions, 
radiation  becomes  a  significant  factor  and  can  not  be  neglected.  Future 
work  to  extend  such  developments  as  Shih's  two-dimensional  radiative 
flux  model  [94,95]  to  the  axisymmetric  case  is  needed  in  order,  to  faci¬ 
litate  meaningful  calculations  for  these  conditions.  However,  before 
such  modifications  to  the  code  are  attempted,  it  should  first  be  used  to 
in  the  analysis  of  heat  transfer  in  engines  configured  with  weak  piston 
barriers  and  using  hypergolic  combustion.  As  the  code  finds  use  in  con¬ 
junction  with  experimental  work,  it  should  become  possible  to  calibrate 
the  code  to  give  better  results.  At  the  same  time,  it  is  expected  that 
the  code  wili  serve  to  help  give  direction  to  follow  on  experimental 
research  in  heat  barrier  piston  engines  and  in  other  low  heat  rejection 
engine  configurations.  When  used  hand  in  hand  with  an  experimental 
effort  it  anticipated  that  the  code  will  serve  to  reduce  time  spent  in 
analysis  while  helping  to  accelerate  the  rate  of  progress  in  these 
important  areas  of  study. 
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ENERGY  EQUATION 


A  general  expression  for  the  energy  equation  is  given  by  [91 
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where 


Q  =  internal  heat  generation, 

Qf  =  radiation  heat  flux  vector. 


4>  =  viscous  dissipation  function, 


For  no  radiation  or  internal  heat  generation  this  reduces  to: 
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which  Ramos  modeled  for  turbulent  flow  [66,69]  as  - 


(A-3) 
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by  assuming  that 


h  =  TTpd)  T. 


In  this  formulation: 


w’'  . 


Equation  A-2  was  used  by  Ramos^^  using  a  22  x  30  mesh  to  give 
results  with  fair  agreement  to  recent  experimental  data  for  two  stroke 
calculations  in  a  motored  IC  engine  with  a  valve  arrangement  very  similar 
to  that  used  in  the  present  analysis. 

However,  the  present  analysis  is  for  a  combusting  engine,  not  a 
motored  one.  Most  researchers  at  this  point,  shift  to  the  use  of  an  energy 
equation  which  is  derived  in  terms  of  total  enthalpy  instead  of  tem¬ 
perature.  However,  because  the  ultimate  objective  of  this  analysis  is  a 
study  of  heat  transfer,  and  for  greater  ease  in  keeping  track  of  tem¬ 
perature  throughtout  the  solution  process,  it  was  decided  best  to  develop 
an  expression  for  the  energy  equation  in  terms  of  temperature.  To  facili¬ 
tate  this  it  was  assumed  that: 

_  n  0 
h  =  Cp  T  +  mi  hi  , 
i  =  l 

where 


Tp  =  5]  mi  Cp^  (T)  , 
i=l 

mi  =  mass  fraction  of  species  i, 

0 

hi  =  heat  of  formulation  of  species  i, 

Cp^  (T)  =  the  specific  heat  at  constant  pressure  of  species  i 
at  temperature  T. 

Substitution  of  h  into  equation  A-2  yields: 


0  -  D  ^  0  0 

--  (p  Cp  T)  +  —  (p  1  mi  hi  )  =  —  (P)  +  ♦  +  V  .(K  V.T) 
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Note  that  the  equation  for  the  conservation  of  an  individual  species 
is  given  by  [92] 


—  (p  mi)  =  7  ‘(pD  7*  mi)  +  p  mi  , 


(A-5) 


where 


mi  =  rate  of  formulation  of  species  i  per  unit  volume. 


Note  that  we  can  multiply  each  expression  in  equation  (A-5)  by  hi 

because  the  heat  of  formation  for  species  i  is  a  scalar  constant.  Summing 

n  .  0 

the  sets  of  resulting  equations  and  solving  for  ^  mi  hi  we  obtain 

i  =  l 

n^og^o  n.o 

p  I  f^i  t^i  =  —  (p  I  f^i  hi)  -  7  •  (p  D  7  •  mi  hi) 
i=l  Dt  i=l  i=l 


With  the  Prandtl  and  Schmidt  numbers  equal  to  one  (and  thus  a  Lewis 
number  of  one),  the  above  expression  can  be  substituted  into  equation 
(A-4)  with  the  following  result: 

(pT)  =  L  [L  (p)  +  *  +  V  .  (K  7  .  T)  -  p  ^  mi  h°  ]  .  (A-6) 

Dt  Cp  Dt  i=l 


It  follows  for  the  case  of  turbulent  reacting  flow  that  the  above 
expression  would  take  the  form  of  equation  (2-6)  with 


.  ,  n.o 

qc  =  -P  mi  hi 

i  =  l 


qr"  =  0 
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APPENDIX  6 


PISTON  DISPLACEMENT  AND  VELOCITY  [65] 


The  equations  for  piston  displacement  Zp  and  piston  velocity  Up  are 
expressed  below  in  terms  of  engine  geometry. 


Zd  =  C  +  4  (l-cos(a)t))  +  L  ^  sin  (uit))^  )  .  (3-1) 


Up  =  uiZ  [sin(a)t)  +  — 


^  s  i  n  ( 2a)t ) 


2  y  l2  -  (t  sin(u»t))2 


(B-2) 


where 


C  =  clearance  (i.e.,  distance  between  TDC  of  the  piston  and  top  of 
cylinder),  (value  of  1.27  cm  used) 


m  =  RPM  iL  (angular  velocity  in  radians  per  second) 
30 

t  =  semi-strqke  length  (value  of  3.81  cm  used) 

» 

L  =  connecting  rod  length  (value  of  20.32  cm  used) 


RPM  =  revolutions  per  minute  (value  of  1900  used) 
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Non-Flow  Region  Temperature  Transformation  Technique 

The  relative  geometric  positions  of  the  near  wall  nodal  points  for 
the  energy  equation  in  the  flow  region  are  invariant  in  relation  to  the 
adjacent  cylinder-piston  wall  nodal  points  in  regions  1  and  along  the 
piston  surface  of  region  3.  As  can  be  seen  from  figure  C-1,  this  is  not 
the  case  for  the  nodal  points  in  the  vicinity  of  the  wall  along  the  side 
of  the  cylinder  of  regions  2  and  3  as  well  as  in  the  area  between  the 
cylinder  and  the  side  of  the  piston.  It  is  desirable  to  be  able  to  con¬ 
tinuously  match  boundary  conditions  between  the  flow  and  non-flow  portions 
of  these  areas  of  the  solution  scheme  even  though  they  are  moving  in  rela¬ 
tion  to  one  another.  It  was  considered  optimum  to  redefine  the  x  posi¬ 
tions  of  the  temperature  nodal  points  in  this  portion  of  region  4  for  each 
new  time  step  so  that  they  always  match  the  x  positions  of  the 
corresponding  flow  field  nodal  points.  If  the  nodal  points  in  the  portion 
of  region  4  corresponding  to  an  i  indicy  value  greater  than  or  equal  to 
izstep  are  allowed  to  move  in  the  x  direction  with  piston  motion,  a  means 
of  accounting  for  the  "old"  values  of  temperature  corresponding  to  the  new 
nodal  points  of  the  new  temperature  control  volumes  is  nece'^sary. 

Adequate  for  the  present  analysis  is  the  simple  linear  technique 
given  below  for  downward  piston  motion: 

If  z(i)  <  zold  (i+1)  - 

.^old  _  -rOldp  ,  .  (z(i)  -  zold(i+l))  +  (zold(i)  -  z(i)) 

>  Pnew  ■  ^ - ^ -  (C-1) 

zold(1)  -  zold(i+l) 

If  z(i)  >  zold(i+l)  - 
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■  zold(i  +  2))  +  T°‘“EEoid  (zo1d(i>l)  -  z(i)) 


old. 


-o1d 


Pnew 


(C-2) 


zold(i+l)  -  zold(i+2) 


where 

=  the  old  temperature  corresponding  to  the  nodal  point 
directly  east  of 

Similarly  for  upward  piston  motion: 

If  z  (i)  >  z  old  (i-1)  - 


^“^"^Pold  -  zold(i-l))  +  (zold(i)  -  z(i)) 

*  Pnpw  “  ~~ 


new 


zold(i)  -  zold(i-l) 


If  z(i)  <  zold(i-l) 


.old 


new 


zold(i-l)  -  zold(i-2) 


'  =  the  old  temperature  corresponding  to  the  nodal  point 


Ww 


old 


old 


directly  west  of  T 


where 


APPENDIX  D 


Solution  Region  Interface  Code  Modifications 


As  can  be  seen  in  Figure  D-1,  the  east  boundary  of  the  control  volume 
for  each  of  the  u-velocity  nodal  points  corresponding  to  the  indicy 
i=izstep  can  be  in  either  region  1  or  region  2  depending  on  the  piston 
position  and  can  actually  move  back  and  forth  between  regions  with  piston 
motion.  Similarly,  the  west  boundaries  of  the  u-velocity  nodal  points 
corresponding  to  i=niz+l  is  also  slaved  to  piston  motion  and  can  move  back 
and  forth  from  region  2  and  region  3.  For  these  two  situations,  the 
values  for  the  expressions  of  x*Up  corresponding  to  the  moving  boundaries 
are  specified  as  follows: 

For  the  east  boundary  of  a  u-velocity  nodal  point  control  volume 
corresponding  to  i=izstep  - 


[.5{zu(izstep+l)  +  zu(izstep))  -  z(istep)] 


(D-1) 


For  the  west  boundary  of  a  u-velocity  nodal  point  control  volume 
corresponding  to  i  =  niz  +  1  - 


[  .5( zu(niz+l)  +  zu(niz))  -  z(istep)] 


(D-2) 


Note  that  for  the  west  boundary  of  the  control  volume  of  the  u- 
velocity  nodal  point  corresponding  to  the  i=izstep+l,  x»up  is  identical  to 
that  given  for  D-1.  Also,  for  the  east  boundary  of  the  control  volume  for 
the  u-velocity  nodal  point  for  i=niz,  D'2  applies.  Obviously,  these  two 
boundaries  also  migrate  between  regions.  Thus,  for  completeness  the  x*up 
values  for  these  boundaries  follow: 


izstep  izstep+1  niz  niz+1 


APPENDIX  E 


RESULTS 

This  Appendix  contains  the  following  Tables  and  Figures: 

Table  E-1  u-velocity  equation  error  analysis 

Table  E-2  v-velocity  equation  error  analysis 

Table  E-3  Turbulent  kinetic  energy  equation  error 
analysi s 

Table  E-4  Turbulent  kinetic  energy  dissipation  rate  equation 
error  analysis 

Table  E-5  Flow  energy  equation  error  analysis 

Table  E-6  Error  analysis  for  combined  u,  v,  k,  c,  T 
equations 

Table  E-7  Error  analysis  for  combined  u,  v  and  k 
equations 

Figure  E-1  Case  3  velocity  vector  plot  for  a  crank  shaft 
angle 

Figure  £-1  Comparison  of  u-velocity  profiles  for  case  3  run 
with  an  identical  run  performed  by  Ramos  and 
Sirignano 

Figure  E-3  Case  3  temperature  profile  for  a  crank  shaft  angle 
of  300° 

Figure  E-4  Case  5  temperature  profile  for  0  =  300°,  run  1 

Figure  E-5  Case  5  temperature  profile  for  0  =  360°,  run  1 

Figure  E-6  Case  5  temperature  profile  for  o  =  372°,  run  1 

Figure  E-7  Case  5  temperature  profile  for  0  =  384°,  run  1 

Figure  E-8  Case  5  temperature  profile  for  0  =  360°,  run  2 

Figure  E-9  Case  5  temperature  profile  for  0  =  366°,  run  2 

Figure  E-10  Case  5  temperature  profile  for  0  =  570°,  run  3 

Figure  E-11  Case  5  velocity  vector  plot  for  0  =  300°,  run  1 

Figure  E-12  Case  5  velocity  vector  plot  for  0  =  360°,  run  1 

Figure  E-13  Case  5  velocity  vector  plot  for  0  =  372°,  run  1 


Figure  t-14  Case  5  velocity  vector  plot  for  0  =  384°,  run  1 

Figure  E-16  Case  5  velocity  vector  plot  for  0  =  360°,  run  2 

Figure  E-16  Case  5  velocity  vector  plot  for  0  =  366°,  run  2 

Figure  E-17  Case  6  velocity  vector  plot  for  0  =  570°,  run  3 

Figure  E-18  Case  5  product  mass  fraction  distribution  for 
0  =  372' 

Figure  E-19  Case  5  specific  heat  values  at  0  =  366°,  run  2 

Figure  E-20  Case  5  fuel  mass  fraction  distribution  for 
0  =  366°,  run  2 

Figure  E-21  Case  5  turbulent  kinetic  energy  distribution  at 
0  =  366°,  run  2 

Figure  E-22  Case  5  turbulent  kinetic  energy  dissipation  rate 
distribution  at  0  =  366°,  run  2 

Figure  E-23  Case  5  effective  viscosity  (ueff)  profile  for 
0  =  366°,  run  2 

Figure  E-24  Case  5  pressure  (P)  versus  crankshaft  angle  (6),  run  2 

Figure  E-25  Case  5  piston  surface  heat  transfer  rate  versus 
crankshaft  angle  (0),  run  2 
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TABLE  E-2  v-velocity  equation  error  analysis 


V  veioclty  tittie=  1.57395 


1 

j 

V  nuffl 

V  ;r.dl 

rnum 

ranai 

3 

3 

-0.13117 

-0,13121 

0.00  301 

0.00304 

3 

4 

-0.13365 

-0.13370 

0.00 306 

0.00310 

3 

5 

-0.  1 3616 

-0. 1 3620 

0.00  j 12 

0.00315 

3 

6 

-0.13360 

-0. 13373 

0.00321 

0.00321 

4 

3 

-0.22701 

-0.22708 

0.00456 

0.00456 

4 

4 

-0,23130 

-0.23138 

0.00459 

0.00464 

4 

3 

-0.23564 

-0.23571 

0.00467 

0.00471 

4 

6 

-0.24002 

-0.24009 

0.00475 

0.00479 

5 

3 

-0.41868 

-0.41881 

0.0081 3 

0.00813 

5 

4 

-0.42660 

-0.42674 

0.00814 

0.00825 

5 

3 

-0.43460 

-0.43474 

0.00803 

0.00836 

5 

6 

-0.44267 

-0.44281 

0.00839 

0.00848 

6 

3 

-0.70617 

-0.70640 

0.01480 

0.01480 

6 

4 

-0.71953 

-0.71977 

0.01499 

0.01499 

6 

5 

-0.73302 

-0.73327 

0.01481 

0.01517 

6 

6 

-0.74664 

-0.74689 

0.01536 

0.01536 

T 

3 

-0,99365 

-0.99400 

0.02304 

0.02  304 

7 

4 

-1.01245 

-1.01281 

0,02331 

0.02331 

7 

5 

-l.03l«3 

-I.O3I8O 

0.02305 

0.02359 

7 

6 

-1.05061 

-1.05097 

0.02366 

0.02388 

8 

3 

-1.08951 

-1.08987 

0.02587 

0.02613 

8 

.« 

-1.11012 

-1.11049 

0.02644 

0.02644 

3 

5 

-1.1309** 

-1.13132 

0.02676 

0.02676 

8 

6 

-1.15196 

-1.15233 

0.02687 

0.02707 

numerical  values  v  velocity 

time: 

1.57895 

4  5  6 

7 

8 

'  -I.u3e-OI  -J.u7e-01  -u.55e-01  -7.68e-01  -1.08e*00  -1.18e*00  -1.19e*00  1.06e*00 

6  -).J5e-02  -1.39e.01  -J.nOe-Oi  -11.436-01  -7.47e-01  -i.05e»00  -1.15e»00  -1.166*00  1.05e*00 

5  -).73«-08  -1.366-01  -2.366-01  -4.356-01  -7.336-01  -1.036*00  -1.13e*00  -1.14e*00  1.046*00 

4  -8.916-02  -1.34e-OI  -2.316-01  -4.27e-0t  -7.206-01  -1.01e*00  -1.116*00  -l.lle*00  1.03e*00 

3  -8.756-02  -1.316-01  -2.276-01  -4.19e-01  -7.06e-01  -9.9“e-01  -1.096*00  -1.096*00  1.02e*00 

2  -3.536-02  -1.296-01  -2.236-01  -4.116-01  -6.93e-01  -9.756-01  -1.07e*00  -1.07e*00  1.016*00 

I  -8.506-02  0.  6*00  0.  6*00  0.  6*00  0.  6*00  0.  6*00  0.  6*00  -1.056*00  1.006*00 

06  a  -.04000  0.  0.10000  0.30000  0.50000  0.90000  1.00000  1.00333 

Ol  a  0.08000  0.12000  0.20767  0.38302  0.64604  0.90905  0.99573  1.00006 


analytic  valu69  v  v6loclty  times  1.57895 


1.2  3  4  5  6  7  8  9  7  = 

J 

7  -9.506-02  -1.436-01  -2.476-01  -4.55e-0I  -7.68e-01  -1.08e*00  -1.18e*00  -1.19e*00  1.06e*00 

6  -9.256-02  -1.396-01  -2.406-0’  -4.436-01  -7.47e-01  -1.05e*00  -1.15e*00  -1.15e»00  1.05e*00 

5  -9.086-02  -1.366-01  -2.366-01  -4.356-01  -7.33e-01  -1.03e*00  -1.13e*00  -1.14e*00  1.04e*00 

a  -8.916-02  -1.346-01  -2.316-01  -4.27e-01  -7.206-01  -1.01e»00  -1.11e*00  -1.11e*00  I.03e*00 

3  -a.756-o;  -1.316-01  -2.27e-01  -4. 19e-01  -7.06e-01  -9.94e-01  -1.09e*00  -1.09e*00  1.02e*00 

2  -8.586-02  -1.296-01  -2.236-01  -4.116-01  -6.93e-01  -9.756-01  -1.07e*00  -1.07e*00  1.01e*00 

I  -8.506-02  0.  6*00  0.  6*00  0.  6*00  0.  6*00  0.  6*00  0.  6*00  -1.06e*00  1.00e*00 

Ol  «  -.04000  0.  0.10000  0.30000  0.60000  0.90000  1.00000  1.00333 

Ol  I  0.08000  0.12000  0.20767  0.38302  0.64604  0.90905  0.99673  1.00005 


TABLE  E-3 


Turbulent  kinetic  enerry  equation  orror  analysis 
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2. 39e*00 

1 .0064.00 

Oi  * 

..C4COO 

0, 

0, 10000 

0. 30000 

0.60000 

0.90000 

1.00000 

1.00333 

Ox  X 

0.C8000 

0. 12000 

0.20767 

0.38302 

0.64604 

0.90906 

0.99673 

1.00006 

analytic  values  turbulence  energy  times  1.5789S 


i 

•  2 

3 

4 

5 

6 

7 

8 

9 

y  = 

J 

7 

1. 15e-01 

1.796-01 

3.32e-Ol 

6.97.-01 

I,39e4.00 

2.25e4.00 

2.58e-400 

2.5964-00 

1.066  4-00 

6 

1. 136-01 

1.75e-01 

3.266-01 

6.86.-01 

1.3764-00 

2.2264-00 

2.5564-00 

2.5664-00 

1. 056*00 

5 

1.1  le-01 

1.72e-01 

3.216-01 

6.79.-01 

1 . 35e*00 

2. 19e4.00 

2.5164-00 

2.5264-00 

1.0464-00 

< 

t.09e-01 

1.696-01 

3. 166-01 

6.6U.-OI 

1.33e4-00 

2. 1664-00 

2.4864-00 

2.4964-00 

1.036  4-00 

3 

1.076-01 

1.666-01 

3. 106-01 

6.611.-01 

1.31e4.00 

2.13e^OO 

2.4464-00 

2.4564-00 

1.0264-00 

2 

i.05e-O1 

1.646-01 

3.056-01 

6.43.-01 

l.29e4-00 

2.  1064-00 

2.4164-00 

2.4264-00 

1.0164-00 

1 

1.03e-O? 

1.616-01 

3.006-01 

6. 33.-01 

1.2764-00 

2,0764-00 

2.  3764-00 

2.396-400 

1.006  4-00 

Oi  * 

-.04000 

0. 

0.10000 

0. 30000 

0.60000 

0.90000 

1.00000 

1.00333 

Ox  t 

0.08000 

0. 12000 

0,20767 

0.38302 

0.64604 

O.90906 

0.99673 

1.00006 
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TABLE  E-4  Turbulent  kinetic  enerry  dissi^^ation 
rate  equation  error  analysis 


«n«rgy  dissipation 


t 

J 

e  null) 

e  anal 

rnun 

ranal 

3 

2 

0.011693 

0.011697 

-.002516 

-.002518 

3 

3 

0.01 1915 

0.01 1919 

-.002592 

-.002589 

3 

4 

0.012140 

0.012144 

-.002662 

-.002662 

3 

5 

0.012366 

0.012J71 

-.002737 

-.002737 

3 

6 

0.012595 

0.012599 

-.002810 

-.002013 

i( 

o 

0.035021 

0.035033 

-.003334 

-.003354 

.* 

4 

3 

0.035687 

0.035699 

-.003434 

-.003428 

U 

4 

0.0^6  359 

0.036  371 

-.003525 

-.003525 

4 

5 

0.  J37037 

0.03^050 

-.003618 

-.003623 

4 

6 

^•.037722 

0.0377  35 

-.0037  32 

-.003725 

«• 

5 

2 

0, 119126 

0. 119166 

-.006478 

-.006478 

5 

3 

u. <21391 

0. 121433 

-.006630 

-.006657 

i 

5 

4 

0. 123678 

0. 123720 

-.006813 

-.006840 

•  -  ■ .  ’ 

5 

3 

0. 125986 

0. 126029 

-.007027 

-.007027 

5 

6 

0. 128315 

0.128359 

-.007219 

-.007219 

6 

2 

0.333906 

0.339023 

-.017955 

-.018040 

6 

3 

0.345350 

0.345470 

-.018625 

-.018533 

6 

4 

0,351355 

0.351978 

-.019028 

-.019028 

6 

3 

0.358421 

0. 358547 

-.019535 

-.019535 

6 

6 

0.365049 

0.365176 

-.0201 14 

-.020054 

7 

2 

0.671031 

0.671267 

-.044494 

-.044494 

fc  - . . , 

7 

i 

0.683790 

0.684033 

-.045672 

-.045672 

7 

4 

0.696670 

0,696918 

-.046075 

-.046875 

7 

5 

0.709672 

0.709924 

-.040106 

-.048106 

‘ 

7 

6 

0.722704 

0.723051 

-.049184 

-.049363 

3 

2 

0.806706 

0.306900 

-.058347 

-.058347 

8 

3 

0.822044 

0.822337 

-.059884 

-.059884 

a 

4 

0,037529 

0.837823 

-.061199 

-.061454 

a 

5 

0.853159 

0.853464 

-.062932 

-.063058 

•  -  *  - 

8 

6 

0.868934 

0.869244 

-.064697 

-.064697 

nutnerlcai  values  energy  dissipation 


L 

1 

J 

7 

S  2 

3 

4 

5 

6 

7 

8 

9 

y  * 

P 

5,7':e-03 

l.28e-02 

3.84e-02 

1.3le-01 

3.72e-01 

7.36e-01 

8.85e-0l 

8.91e-01 

1 .06e-*-00 

6 

5e53e-03 

1.26e-02 

3.77e-02 

1.28e-01 

3.65e-01 

7.23e-0l 

8.69e-01 

8.75e-01 

1.05e*00 

5 

5.50e-03 

t.24e-02 

3.70e-02 

1.26e-01 

3.58e-01 

7. lOe-01 

8.53e-Oi 

8.59e-Ol 

1 .04e«-00 

« 

5e40«-03 

1.2l«-02 

3.64e-02 

1.24e-01 

3.52e-01 

6.97e-01 

8.38e-Ol 

8.43e-01 

1.03e»00 

3 

5.30e-03 

l,l9e-02 

3.57e-02 

1.2le-0t 

3.«5e-01 

6,84e-Ol 

8.22e-0i 

8.28e-01 

1.02e*00 

2 

5.20e-03 

1. 17e-02 

3.5Q«-02 

1.19e-01 

3.39e-01 

6.71e-01 

8.07e-01 

8.12e-01 

1 .01e«-00 

1 

5.l0e-03 

1. l5e-02 

3.44e-02 

I.i7e-0l 

3.33e-Oi 

6.59e-01 

7.92e-01 

7.97e-01 

t.00e*00 

t 

Cl  s 

-.04000 

0. 

0. 10000 

0.30000 

0.60000 

0.90000 

1.00000 

UOO333 

m 

Ox  * 

0.08000 

0. 12000 

0.20767 

0.38302 

0.64604 

0,90906 

0.99673 

1.00006 

■  r 


analytic  values  energy  dissipation  timer  1.S7895 


L 

1 

*  2 

3 

4 

5 

6 

7 

8 

9 

y  = 

1 

7 

5.70e-03 

1.28e-02 

3.84e-02 

1,31e-01 

3.72e-01 

7.36e-01 

8.85e-01 

8.9l«-Ol 

1.06e<^00 

6 

5.60e-03 

1.26e-02 

3.77e-02 

1.28e-01 

3.65e-01 

7.23e-01 

8.69e-01 

0.75«-O1 

1.05e*00 

5 

5.50«-03 

1.24e-02 

3.71«-02 

1.26e-01 

3.59e-01 

7.10e-01 

8.53e-01 

8,59«-01 

1.04e«00 

V 

5.40e-03 

1.21e-02 

3.64e-02 

l.24e-01 

3.52e-01 

6.97e-01 

8.30«-O1 

8.43e-01 

1.03e4-00 

3 

5.)0e-03 

1. 19e-02 

3.57e-02 

1.21C-01 

3.45e-01 

6.84e-0l 

d.22e-01 

6.28e-01 

1.02e4>00 

2 

5.20e-03 

1.17e-02 

3.50e-02 

1. 19e-01 

3.39e-0l 

6.71e-01 

8.07e-01 

8. 12e-01 

1  .Ole^OO 

1 

5.l0e-03 

1.15S-02 

3.44e-02 

1. 17e-01 

3.33e-Ol 

6.59e-Ol 

7.92e-01 

7.97e-01 

1.00e-*-00 

t 

-.04000 

0. 

0. 10000 

0,30000 

0.60000 

0.90000 

1.00000 

1.00333 

t 

o.oaooo 

0. 12000 

0.20767 

0,38302 

0.64604 

0.90906 

0.99673 

1.00006 

'“'-■A'i 

t-'  iv'.Vi 
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TABLE  E-5  Flow  energy  Equation  error  analysis 


temperature  times  1.57395 


t 


■ 


1 

J 

t  nun 

t  anal 

rnum 

ranal 

3 

2 

0.023309 

0.023396 

-.002394 

-.002403 

3 

3 

0,02383« 

0.023841 

-.002499 

-.002503 

3 

4 

0,024283 

0.024290 

-.002607 

-.0Hcb07 

3 

5 

0.024736 

0.024744 

-.002710 

-.002714 

3 

6 

0.025193 

0.025201 

-.002814 

-.002825 

4 

•5 

0.070049 

0.070072 

-.003917 

-.00  3929 

4 

3 

0,071331 

0.071405 

-.004050 

-.004073 

4 

4 

0.072726 

0.072750 

-.004187 

-.004221 

4 

5 

0.074083 

0.074108 

-.004352 

-.004375 

4 

6 

0.075453 

0.075473 

-.004533 

-.0049  33 

5 

2 

0.238276 

0.238357 

-.010099 

-.010211 

5 

3 

0.242307 

0.242889 

-.010474 

-.010529 

5 

4 

0.247 381 

0.247465 

-.010800 

-.010855 

5 

5 

0.251998 

0.252083 

-.01 1 1 36 

-.01 1 190 

5 

6 

0.256657 

0.256744 

-.011426 

-.011533 

6 

2 

0.677083 

0.678114 

-.031359 

-.031545 

6 

3 

0.690772 

0.691010 

-.032424 

-.032424 

6 

4 

0.703734 

0.704027 

-.033505 

-.033322 

6 

5 

0.716919 

0.717166 

-.033878 

-.034240 

6 

6 

0.730174 

0,730426 

-.034939 

-.035173 

7 

2 

1.342205 

1.342668 

-.075901 

-.075622 

7 

3 

1.36772$ 

1 . 368202 

-.077311 

-.077587 

7 

4 

1.393690 

1.393976 

-.079044 

-.079592 

7 

5 

1.111969$ 

1.419991 

-.081637 

-.081637 

7 

6 

1.ii6$762 

1.446246 

-.083006 

-.083722 

3 

2 

I. 6 13586 

1.614141 

-.097298 

-.007557 

8 

3 

1.66II26II 

1.644838 

-. 100298 

100041 

8 

4 

1.675239 

1.675823 

.. 102317 

-.102572 

8 

5 

1.706502 

1.707098 

-.104900 

-.105153 

8 

6 

1.738055 

1.738661 

-.107450 

-. 107783 

numerical  values  temperature  time*  1,57095 


1 

s  2 

3 

4 

5 

6 

7 

8 

9 

y  s 

1. 14e-02 

2.57e-02 

7,69e-02 

2,61e-01 

7.44e-0l 

l.47e+00 

1.7764.00 

1.7864-00 

1.0664-00 

6 

1.  l2e-02 

2.52e-02 

7.55e-a2 

2.57e-0t 

r.30e-0l 

I.45e-«>00 

1.74e4.00 

1.7564-00 

1.0564-00 

5 

1. iCe-02 

2,«7e-02 

7.uie-02 

2.52e-Ol 

7.l7e-01 

1 .42e^00 

1.7le^00 

1.7264-00 

1.0464-00 

S 

l.Ode-02 

2.43e-02 

7,27e-02 

2.47e-01 

7.04e-01 

1.39e4.00 

1 ,68e*00 

1.6964-00 

1.0364-00 

3 

1.06e-02 

2,38e-02 

7.14e-02 

2.43e-01 

6.91e-01 

1. 3704.00 

1 .64e4.00 

1.6664-00 

1.0264-00 

2 

l.oae-02 

2.34e-02 

7.00e-02 

2,38e-01 

6.78e-0l 

1.34e4.00 

1 .6 1e*00 

1.6264-00 

1.0164-00 

1 

1.02e-02 

2.30e-02 

6.88e-02 

2.34e-01 

6.65e-01 

1. 3264.00 

1.58e+00 

1.5964-00 

1.006400 

Ox  s 

-.04000 

0. 

0. 10000 

0.30000 

0.60000 

0,90000 

1.00000 

UOO333 

0*  s 

0.08000 

0.12000 

0.20767 

0.38302 

0.64604 

0.90906 

0.99673 

1.00006 

analytic  values  temperature 


times  1,57895 


1 

J 

7 

S  2 

3 

4 

5 

6 

7 

8 

9 

y  « 

f 

1. 14e-02 

2.57e-02 

7.69e-02 

2.61e-Ol 

7.«4e-0I 

1.476400 

1.776400 

1,786400 

1.006400 

6 

1.12e-02 

2.52e-02 

7.556-02 

2.572-01 

7.306-01 

1.456400 

1.746400 

1.756400 

1.056400 

5 

1.10e-02 

2.47e-02 

7,41e_02 

2-52e-01 

7.176-01 

1.426400 

1.716400 

1,726400 

1.046400 

«  ■ 

4 

1.086-02 

2.43e-02 

7.276-02 

2.47e-01 

7.04e-01 

1.396400 

1.606400 

1.696400 

1.036400 

3 

1.06e-02 

2.38e-02 

7. 1«e-02 

2.436-01 

6.91e-01 

1.376400 

1.646400 

1.666400 

1.026400 

2 

1.04e-02 

2.346-02 

7.016-02 

2.386-01 

6.78e-01 

1.346400 

1.616400 

1.626400 

1.016400 

1 

U02e-02 

2.30e-02 

6.88e-02 

2.34e-01 

6.65e-01 

1.326400 

1.586400 

1.596400 

1.006400 

Ox  s 

-.04000 

0. 

0. 10000 

0.30000 

0.60000 

0.90000 

1.00000 

1.00333 

Oz  s 

0.08000 

0. 12000 

0.20767 

0,38302 

0.64604 

0.90906 

0.99673 

1.00006 
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TABLE  E-6  Error  analysis  for  combined 
u,  V,  k,  e,  T  equations 


tlM*  1.57895 


tlM«  1.57895 


tiirtHilcnc*  •ntrtir 


■"*  ■  ‘’.V* 


j 

u  nun 

<J  4n4l 

rnu« 

ranjl 

l 

i 

V  ftun 

V  «n4l 

rnun 

ranjit 

1 

3 

k  nun 

k  an«l 

rnun  ranal 

? 

n.04.'i79 

0.04,?993 

-.001225 

-.001225 

3 

3 

-0. 1 31 17 

-0.13'.'' 

0.00  30 ' 

0.00  104 

3 

2 

0. 16)471 

0. 163S26 

-.000841  -.008041 

3 

0.04 

0.0434110 

-.QQ1220 

-.0012)5 

3 

ii 

-Q. 1)365 

-0.1)370 

0.00 )01 

0.00)10 

3 

3 

0. 166414 

0. '66470 

-.00906?  -.C09062 

« 

0.043792 

0.043007 

-.001246 

-.001246 

3 

5 

-0. 1 J6t6 

-0.13620 

0. 00315 

O.OOJI5 

3 

4 

0. 16930) 

0. 169441 

-.009207  -.009207 

■  - 

5 

0.344199 

9.344214 

-.001249 

-.001256 

3 

6 

-0. '3069 

-0.13873 

0.00)16 

0.00)21 

3 

5 

0.172)79 

0.1724)8 

-.009518  -.009518 

6 

0.044606 

0.044621 

-.001257 

-.001266 

4 

3 

-0.2270' 

-0.22708 

0.00447 

0.00456 

3 

6 

0.175401 

0.  175461 

-.009048  -.CO9755 

2 

0.1 J0667 

0.1 39714 

-.002150 

-.002150 

4 

4 

-0.23130 

-0.23'38 

0.00459 

0. 00464 

4 

2 

0. 304770 

0.304071 

-.011001  -.011047 

3 

0. 140989 

0. 141036 

-.002176 

-.002176 

4 

5 

-0.23564 

-0.2357' 

0.00467 

0.00471 

4 

3 

0. 310070 

0. 31017) 

-.012003  -.0121)0 

It 

0. 14?3i  I 

0. 142)50 

-.002203 

-.00220) 

4 

6 

-0.24002 

-0.24009 

0.00471 

0.00479 

4 

4 

0.315416 

0.315521 

-.012418  -.012418 

5 

0.1436)3 

0. 143600 

-.002162 

-.002230 

5 

3 

-0.41066 

-0.41801 

0.0U802 

0.0001 3 

4 

5 

0. 320000 

0.320915 

-.01262)  -,01271) 

6 

0. 144954 

0. 14S003 

-.002257 

-.002257 

5 

4 

-0.42660 

-0.42674 

0.00014 

0.00825 

4 

6 

0. J26246 

0.326354 

-.013075  -.013015 

2 

0.423089 

0.4240)0 

-.003822 

-.00)022 

5 

5 

-0.4)460 

-0.43474 

o.f'oaiu 

O.OO0 36 

5 

2 

0.642750 

0.64^96? 

-.021524  -.021524 

3 

0. 427900 

0.42004) 

-.00 3863 

-.00  3863 

5 

6 

-0.44267 

-0.44261 

0.00839 

0.006  <10 

5 

3 

0.65 3288 

0.65)506 

-.021905  -.022016 

II 

0.431911 

0.4)2055 

-.00)705 

-.00  1904 

6 

3 

-0.70617 

-0.70640 

0.0»«43 

0.01400 

5 

4 

0.663912 

0.6641)3 

-.022408  -.022517 

5 

0.4)5923 

0.416060 

-.00)749 

-.00)946 

6 

0 

-0.71953 

-0.7«977 

0.01499 

0.01499 

5 

5 

0.674619 

0.674044 

-.02)1)8  -.02)029 

. 

6 

0. 4  399)6 

0.440001 

-.003923 

-.00)900 

6 

5 

-0.73302 

-0.73327 

0.01499 

0.01517 

5 

6 

0,605411 

0.685638 

-.023552  -.02)552 

2 

0.960029 

0.960)47 

-.005062 

-.005509 

6 

6 

-0.74664 

-0.74689 

0.01536 

0.915  36 

6 

2 

1.288175 

1.20060) 

-.048640  -.040865 

•  * 

3 

3  977107 

0. 977512 

-.005006 

-.005528 

7 

3 

-0.99365 

-0.99400 

O.OT.’OO 

0.0? ’04 

6 

3 

1.  )0799l 

1.300319 

-.049959  -.049959 

« 

0  )»6»u0 

0.986676 

-.005109 

-.005540 

7 

4 

-'.0'245 

-t.0'2|l 

0.02  3"5 

0.0233* 

6 

4 

1.  3277  30 

1.)?e\76 

-.051073  -.05107) 

•* 

5 

0.995500 

0.995040 

-.005569 

-.005568 

7 

5 

-'.0)14) 

-1.0J'80 

0.02 305 

0.02 159 

6 

5 

1. 347721 

1.340173 

-.052390  -.052209 

6 

1.004673 

1.005004 

-.005506 

-.005500 

7 

6 

-1.0506' 

-1.05097 

0.02388 

0.02308 

6 

6 

1. 367862 

1 . 368)12 

-.053)68  -.053)60 

2 

1.453068 

1.454)40 

-.005607 

-.006127 

a 

J 

-'.0095' 

.'.00987 

0.0261) 

0.0261) 

7 

2 

2.099759 

2. 100449 

-.099271  -.099029 

w  - 

3 

1.467612 

1.460104 

-.005660 

-.006096 

0 

0 

-'."0'2 

.'.It049 

0.02619 

0.02644 

7 

3 

2. 1)0200 

2.130909 

-.102570  -.102016 

« 

1.401372 

1.401067 

-.005848 

-.006064 

0 

5 

•'.'3094 

-1.i3'J2 

0.02574 

0.02676 

7 

4 

2.  16085' 

7. 161568 

-.103697  -.10«2«5 

,* 

5 

'.495132 

1.4956)1 

-.005819 

-.00603) 

8 

6 

-'.'5196 

-1.15233 

0.0268T 

0.02707 

7 

S 

2. 191700 

2.192426 

-.106516  -.106516 

•  ' 

6 

1.500096 

1.509)94 

-.005437 

-.006001 

7 

6 

2.222741 

2.22)402 

•.1O9108  -.100030 

8 

2 

2.407218 

2.407998 

-.12)704  -.120223 

8 

) 

2.411588 

2. ••2390 

-.126405  -.'26919 

8 

0 

2.076178 

2.477000 

-.1)0176  -.129656 

8 

5 

2.510998 

2.511827 

-.13»960  -.132»65 

_ 8 

6 

2.5460)3 

2.546872 

-.135116  -.135316 

. 

«ntrc)r  4imp«ti«n  '.57895 


llm*  '.57895 


1 

J 

«  nw 

•  •nal 

rnufl 

rOMt 

1 

J 

t  nu* 

t  6ft*l 

rnua  ranal 

3 

2 

0.011693 

O.01169T 

-.002516 

-.002510 

) 

2 

0.02)389 

0.023396 

-.00240)  -.00240) 

c. 

3 

3 

0,011915 

0.011919 

-.002592 

-.002509 

) 

) 

0.0238)4 

0.02)041 

-.00250)  -.00250) 

s*. 

3 

0 

0.012140 

O.O121«0 

-.002660 

-.002662 

) 

0 

0.024283 

0.024290 

-.002598  -.002607 

3 

5 

0.012)66 

0.012)73 

-.0027)7 

-.002757 

5 

5 

0.0247)6 

0.024744 

-.002714  -.002714 

3 

6 

0.012595 

0.012599 

-.00281) 

-.00201) 

) 

6 

0.025193 

0.025201 

-.002819  -.002825 

4 

2 

0.035021 

0.035033 

-.00)328 

-.00)))4 

0 

2 

0,070049 

0.070072 

-.00)894  -.00)929 

m 

0 

3 

0.035687 

0.035699 

-.00)422 

..00)020 

0 

5 

0.071)81 

0.071005 

-.0040)8  -.004073 

4 

4 

0.036359 

0.036)71 

-.00)519 

-.00)525 

4 

0 

0.072726 

0.072750 

-.004107  -.004221 

4 

5 

0.037038 

0.0)7050 

-.00)629 

-.00)62) 

4 

5 

0.074083 

0.074108 

-.004)64  -.004)75 

4 

6 

0.0)7723 

0.037735 

-.003725 

-.00)725 

0 

6 

0.07545) 

0.075478 

-.004503  -.004533 

5 

2 

0.119127 

O.II9166 

-.006451 

-.006070 

5 

2 

0.2)8276 

0.238357 

-.010267  -.010211 

5 

3 

0.121392 

0.1214)) 

-.006657 

..006657 

5 

) 

0.242807 

0.242009 

-.010474  -.010529 

5 

0 

0. 12)679 

0.12)720 

-.00683) 

. .006840 

5 

« 

0.247381 

0.247465 

-.010800  -.010855 

.  ^ 

5 

5 

0.125987 

0.126029 

-.007027 

..007027 

5 

5 

0.25 '998 

0.252083 

-.011190  -.011190 

.  * 

5 

6 

0. 128)16 

0. 128359 

-.007219 

-.007219 

5 

6 

0.256657 

0.256740 

-.011498  -.0115)3 

«* 

6 

2 

0.338911 

0.3)9023 

-.018094 

-.0180*0 

6 

2 

0.677882 

0.678110 

-.031545  -.031545 

6 

3 

0.345355 

0. 345470 

-.01853) 

-.0185)) 

6 

3 

0.690772 

0.691010 

-.032424  -.032424 

*  - 

6 

0 

0.351861 

0.351978 

-.0189)7 

-.019020 

6 

0 

0.703704 

0.700027 

-.0)3322  -.03)322 

6 

5 

0.350427 

0.358547 

-.019535 

-.0195)5 

6 

5 

0.716919 

0.717166 

-.03)878  -.0)4200 

6 

6 

0.365055 

0.365176 

-.020054 

-.020054 

6 

6 

0.7)0174 

0.7)0026 

-.035058  -.035178 

7 

2 

0.671045 

0.671267 

-.044)50 

.,04*094 

7 

2 

1.342205 

1.342668 

-.075)42  -.075622 

7 

3 

0.683805 

0.684033 

-.0455)) 

-.045672 

7 

) 

1.367725 

1.368202 

-.077587  -.077587 

", 

7 

0 

0.696686 

0.6969*8 

-.046875 

-.046815 

7 

0 

1. 393490 

1.393976 

-.079044  -.079592 

7 

5 

0.709688 

0.709920 

-.047970 

-.040106 

7 

5 

1.419495 

1.419991 

-.080822  -.08)637 

‘  *- 

7 

6 

0.722810 

0.723051 

-.049542 

-.049)6) 

7 

6 

1.445741 

1.4*6246 

-.00)722  -.00)722 

8 

2 

0.806726 

0.006990 

-.050477 

-.050)47 

0 

2 

1.613586 

1.614141 

-.097290  -.097557 

8 

3 

0.822065 

0.822)37 

-.059755 

-.059004 

8 

) 

1.640265 

1.6400)8 

-.100041  -.100001 

8 

0 

0.8)7550 

0.837828 

-.061450 

-.061454 

0 

4 

1.675239 

1.67582) 

-.102)17  -.102572 

*  •. 

8 

5 

0.85J181 

0.05)464 

..06)051 

-.06)050 

0 

5 

t. 706502 

1.707098 

-.100900  -.10515) 

8 

6 

0.868957 

0.869244 

-.0605)1 

-.064697 

0 

6 

1.7)8056 

1.7)8661 

-.107117  -.10770) 

-  f 


ro  ^ 
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TABLE  E-6  Error  analysis  for  corbined 

u,  V,  k.  T  equations  (cont.) 


nunerical  values  tirnes  1.3739*? 


u  velocity 


1  s  3 

u 

5 

6 

7 

8 

9 

y  = 

J 

l.o-e-oa 

4.30e-02 

1 .46e-01 

4.44e-01 

1.0le-*-00 

1.520-00 

1 .68e*00 

1.066400 

a  ?.o6«-oa 

4.-46e-02 

t.45«-0l 

4.40e-01 

I.OOe^OO 

1.3Je^00 

1.66e-00 

1.056400 

3  l.C3e-.03 

4.42e-02 

1 .4Ue-01 

4. 36e-Ol 

9.96e-01 

1.50e-00 

1 . 65e-00 

1.046-00 

4  i.C<4e«U3 

4. 38e-02 

1  .U2e-01 

4.32e-01 

9.86e-01 

1 .48e-00 

1 .63e-00 

1 .03e-00 

3  i.03e-02 

4. 34e-02 

1.41e-01 

4.28e-0l 

9.77e-01 

1.47e-00 

1.62e-00 

1.026-00 

2  t.03e-'^2 

4.30e-02 

1 .40e-01 

4.24e-01 

9.68e-01 

1.45e+00 

1 .60e-00 

1 .016-00 

1  l.Q2e-02 

4,26e-02 

1.38e-01 

4.20e-01 

9.590-01 

1 . 44e-00 

1 .59e-00 

1  .OOe-00 

Oi  s  -.04000 

0.03000 

0.20000 

0.43000 

0.75000 

0.95000 

1.00333 

oz  s  o.oaooo 

0.16384 

0.29335 

0.31453 

0.77755 

0.95289 

1.00006 

V 

velocity 

1  »  2 

J 

7  -9.50«-Q2 

3 

4 

5 

6 

7 

a 

9 

y  = 

-l.43e-Ol 

-2.47e-01 

•4.35e-0l 

-7.68e-01 

-1 .08e-00 

-1.18e-00 

-1 . 19e-00 

1 .06e-00 

6  -9.25«-02 

-l.39e-01 

-2.40e-01 

-4.43e-01 

-7.47e-01 

-1.05e>00 

-J.  15e-00 

-).  l6e-00 

1.05e-00 

S  -9.08«-02 

-l.36e-Ot 

-2.36e-01 

-4. 33e-01 

-7.330-01 

-1,03e-00 

-1 . 13e-00 

-1 . 14e-00 

1 .046400 

4  -a.91e-02 

-1 . 34e-Oi 

-2.31«-01 

-4.27e-0l 

-7.20e-01 

-1  .Ole-00 

-1 . 1  le-00 

-1 .  lle-00 

1 .03e-00 

3  -a.75«-’:2 

-l.3i«-0i 

-2.27e-0t 

-4. 19e-01 

-7.06e-0I 

-9.94e-01 

-!  .09e-00 

-1 .09e-00 

1 .02e-00 

2  -3.38e-02 

-t.29«-0l 

-2.23e-01 

-4. 1 te-01 

-6.93e-01 

-9.75e-01 

-1 .07e*00 

-1 .07e-00 

1 .01e-00 

J  -a.30e-02 

0.  e*00 

0.  e^OO 

0.  e*00 

0.  e^OO 

0.  e-00 

0.  e-00 

-1 .06e-00 

1.006400 

't  a  -.04000 

Q. 

0.  toooo 

0. 30000 

0.60000 

0.90000 

1.00000 

1.00333 

vZ  a  0,03000 

0.  12000 

0.20767 

0.38302 

0.64604 

0.90906 

0.99673 

1.00006 

turbulence  energy 

t  s  2 

J 

7  M5«-01 

3 

4 

5 

6 

7 

8 

9 

y  s 

1.79e-01 

3,32e-01 

6,97e-01 

1.39e*00 

2.25e*00 

2.3de-00 

2,59e-00 

1.066400 

6  I.t3#-01 

3,26e-0} 

6.a5«-o? 

1,37e*00 

2.22e-00 

2.55e4.00 

2.566400 

1.056400 

3  i,ne-0i 

1.72«-01 

3.2t«-01 

6.73«-01 

1.35€^00 

2.19e-00 

2.31e*00 

2.526400 

1.046400 

u  i.09«<-ot 

1 .69«-01 

3«15«-01 

6.64e-0l 

1.33e*00 

2.16e4-00 

2.486400 

2.496400 

1.036400 

3  l,07e-0t 

f,36«-Of 

3.  !0e-0i 

6.53«-0t 

t.31e*00 

2. I3e400 

2.446400 

2.436400 

1.026400 

2  1,03«-01 

i.63e-Ol 

3.03e-01 

6.43e-01 

1.29e^00 

2.10e-00 

2. 418400 

2.426400 

1.016400 

1  l.03e-Oi 

l.6le-01 

3.00e-01 

6.33«-Ol 

1.27e*00 

2,07e-00 

2.376400 

2.396400 

1. 006400 

Ot  a  -.QUOOO 

0. 

0. 10000 

0.30000 

0.60000 

0.90000 

1,00000 

1.00333 

Os  S  0,08000 

0, 12000 

0,20767 

0.38302 

0.64604 

0.90906 

0.99673 

1.00006 

en6rgy  dissipation 

1 

a  2 

3 

4 

5 

6 

7 

8 

9 

y  5 

j 

7 

5.706-03 

1.286-02 

3.846-02 

1,316-01 

3.726-01 

7,36e-01 

3.856-01 

8.91e-01 

1.06e-*-00 

6 

5.606-03 

1.266-02 

3.776-02 

1.286-01 

3.656-01 

7.236-01 

8.69e-01 

8.75e-01 

1.05e<»>00 

5.506-03 

1.24«.02 

3,706-02 

1.266-01 

3.506-01 

7,10e-01 

8.53e-0l 

3.596-01 

1 .04e4>00 

« 

5.406-03 

1.216-02 

3.646-02 

1.246-01 

3.526-01 

6.97e-01 

8.38e-0l 

8.43e-01 

1.036400 

3 

5.306-03 

1.196-02 

3,576-02 

1.216-01 

3.456-01 

6.846-01 

8.22e-Ol 

8.28e-01 

1.026400 

2 

5.206-03 

1.176-02 

3.506-02 

1.196-01 

3.396-01 

6,71e-01 

8.07€-01 

0.126-01 

1.016400 

1 

5.106-03 

1.156-02 

3,44e-02 

1,176-01 

3.336-01 

6.59e-01 

7.92e-0l 

7.976-01 

1.006  4-00 

Os  a 

-.04000 

0. 

0.10000 

0.30000 

0.60000 

0.90000 

1.00000 

1.00333 

Os  a 

0.08000 

0.  12000 

0.20767 

0.38302 

0.64604 

0.90906 

0.99673 

1.00006 

t6mp6rature 

1 

a  2 

3 

4 

5 

6 

7 

8 

9 

y  * 

J 

7 

f.  146-02 

2.576-02 

7.696-02 

2.616-01 

7.44e-0J 

1.476400 

1.776400 

1.786400 

1.066400 

6 

1.126-02 

2.526-02 

7.556-02 

2,576-01 

7.30e-01 

1.456400 

1.746400 

1.756400 

1.056400 

5 

1.106-02 

2.476-02 

7,416-02 

2.526-01 

7.176-01 

1.426400 

1.716400 

1.726400 

1.046400 

4 

1.086-02 

2.436-02 

7.276-02 

2.476-01 

7.046-01 

1.396400 

1.686400 

1.696400 

1.036400 

3 

1.066-02 

2.386-02 

7. 1«e-02 

2.436-01 

6.916-01 

1.376400 

1. 64-400 

1 .666400 

1.026400 

2 

1.046-02 

2.346-02 

7.006-02 

2.386-01 

6,78e-01 

1.346400 

1.616400 

1.626400 

1.016400 

1 

1.026-02 

2.306-02 

6.886-02 

2.346-01 

6.656-01 

1.326400 

1.586400 

1.596400 

1.006400 

“t  a 

-.04000 

0. 

0.10000 

0.30000 

0.60000 

0.90000 

1.00000 

1.00333 

'Z  a 

0.08000 

0.12000 

0.20767 

0.38302 

0,64604 

0.90906 

0.99673 

1.00006 

TABLE  E-b 


Error  analysis  for  combined 
u,  V.  k,  ‘I  T  equations  (cont.) 


analytic  values 


tlme=  1.57895 


u  velocity 


L  -  3 

J 

’  l.07e-02 

4 

5 

6 

7 

8 

9 

y  = 

4.50e-Q2 

1.46e-01 

4.44e-01 

i.Ole^OO 

1 .52e*00 

1.686*00 

1 .06e*00 

6  1.06e-02 

4.46e-02 

l,45e-01 

4.40e-01 

1 .01e«^00 

1.5164.00 

1.666*00 

1 .056*00 

5  l.a5«-02 

4.42e-02 

1 .4Ue*01 

4. 36e-01 

9.96e-01 

1.50e+00 

1.656*00 

1 .04e*00 

4  l.Q4«-02 

4.i8«-Q2 

1.42e-01 

4.32e-0\ 

9.87e-01 

1.486  4.00 

1 .636*00 

1 .036*00 

3  i.03e-02 

4. 34e-02 

1.41e-01 

4.23e-01 

9.78e-01 

1 .47e*00 

1.626*00 

1 .026*00 

2  1.03e-02 

4. 30e-02 

1.40e-01 

4.24e-01 

9.68e-Ol 

1.45e*00 

1 .60o*00 

1.016*00 

1  l.02e-02 

4.26e-02 

1.38e-01 

4.20e-0t 

9.59e-01 

1.44e*00 

1.596*00 

1.00e*00 

a  -.04000 

0.05000 

0.20000 

0.45000 

0.75000 

0.95000 

1.00333 

a  0.08000 

0.  16384 

0.29535 

0.51453 

0.77755 

0.95289 

1.00006 

V 

velocity 

1  a  2 

J 

7  -9,50«-02 

3 

4 

5 

6 

7 

8 

9 

y  = 

-t.43e-ai 

-2.47e-0I 

-4.55e-01 

-7.68e-01 

-1.0864.00 

-1 . 18e*00 

-1 . 196*00 

1.066*00 

6  -9.25e-02 

-t.39«-01 

-2.40e-01 

-4.43e-Ot 

-7.47e-01 

-1.056*00 

-1 . 156*00 

-1.16e*00 

1.05e*00 

5  -9.08e-02 

-l.36e-Ol 

-2.36e-Ol 

-4.35e-01 

-7.33e-01 

-1.03e*00 

-1.13e*00 

-1 . 14e*00 

1.04e*00 

4  .d.9le-02 

-1.34e-01 

-2.31e-01 

-4.27e-01 

-7.20e-01 

-1.016*00 

-1.116*00 

-1 . 1 1e*00 

1.03e*00 

3  -a.75e-02 

-i.3t«-0i 

-2.27e-01 

-4.l9e-01 

-7.06e-01 

-9.946-01 

-1 .096*00 

-1.096*00 

1.02e*00 

2  -8.53e-02 

-l.29«-01 

-2.23e-01 

-4. lle-Ol 

-6.93e-01 

-9.756-01 

-1.07e*00 

-1.076*00 

T.01e*00 

1  -3.50e-02 

0.  e^OO 

0.  e»00 

0.  eoOO 

0.  e>00 

0.  6*00 

0.  e*00 

-1  ,Loe*00 

1.006*00 

Oi 

a  *.04000 

0. 

0. 10000 

0.30000 

0.60000 

0.90000 

1.00000 

K00333 

Oz 

a  0.08000 

0. 12000 

0.20767 

0.38302 

0.64604 

0.90906 

0.99673 

1.00006 

turbulence  energy 

1  a  2 

J 

7 

3 

4 

5 

6 

7 

8 

9 

y  s 

1.79«-01 

3.32e.01 

6.97e-01 

1.39e»00 

2.25e*00 

2.58e*00 

2.596*00 

1.06e*00 

6  1.13«-01 

1.75«-01 

3.26«-01 

6.86e-0l 

1.37e*00 

2.22e*00 

2.55e*00 

2.56e*00 

1.056*00 

5  MU-01 

1.72e-Ot 

3.2U-Ot 

6.75e.0l 

l.35e*00 

2. 19e*00 

2.5le*00 

2,52e*00 

1.04e*00 

4  1.09«-01 

t.69«-01 

3.16e-Ot 

6.64e-01 

1.33e^00 

2.16e*00 

2.U8e*00 

2.49e*00 

l.03e*00 

3  l.07t-01 

1.66«-01 

3. tOe-01 

6.5«e-Ot 

1.3U*00 

2.136*00 

2.44e*00 

2.45e*00 

1.026*00 

2  t-05«-0t 

1.64«*01 

3.05«-01 

6,43e-01 

1.29e*00 

2. 10e*00 

2.4le*00 

2.42e*00 

1.01e*00 

1  l.03«-0l 

1.6U-01 

3.00e-01 

6.33e-Ol 

1.27C4.00 

2.076*00 

2.376*00 

2.39e*00 

1.00e*00 

T* 

a  -.04000 

0. 

0. 10000 

0.30000 

0.60000 

0.90000 

1.00000 

1.00333 

'Z 

a  0.08000 

0.  12000 

0.20767 

0.30302 

0.64604 

0.90906 

0.99673 

1.00006 

energy  dissipation 

1 

J 

7 

•  2 

3 

Q 

5 

6 

7 

8 

9 

y  = 

5.70*-03 

1.23e-02 

3.84e-02 

1.3U-01 

3.72e-01 

7.36e-01 

8.056-01 

8.916-01 

1.06e*00 

6 

5.60e-03 

1.26e-02 

3.776-02 

l.28e-01 

3.65e-Ol 

7.236-01 

8.696-01 

8.756-01 

1.056*00 

5 

5.50e-0j 

i,24e-02 

3.7U-02 

1.26e-0l 

3.59e-01 

7.  lOe-01 

8.536-01 

8.59e-01 

1.04e*00 

e 

5.40e-03 

1.2te-02 

3.64e.02 

1.24e-oi 

3.52e-01 
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Error  analysis  for  combined  u,  v  and  k  equations  (cont.) 
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Figure  E-25  Case  5  piston  surface  heat  transfer  rate  versus  crankshaft 
angle  (e),  run  2  (rate  through  1  radian  of  the  piston's 
surface  area) 
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